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PREFACE 


NASA Lewis Research Center developed and has been conducting research on 
an enlarged engine structures program since 1979. Development of advanced 
methods of nonlinear structural analysis of engine components is a significant 
part this enlarged engine structures program and is an integrated research 
effort involving Lewis, industry, and the university community. 

A two-day workshop was held at Lewis Research Center on April 19 and 20, 
1983, to report recent progress in nonlinear structural analysis for engine 
structures. The workshop was organized into three sessions as follows: 

Session I - New Concepts/Formulations, Session II - Algorithms/Convergence, 
and Session III - Inelastic Analysis and Interactive Elements. 

New concepts/formulations include (1) the slave finite-element formula- 
tion for space and time where interpolation polynomials (amenable to explicit 
integration) are used to express all variables entering the formulation over 
the element domain, (2) new variational principles leading to new formulations 
for hybrid stress finite elements that possess ideal characteristics, such as 
minimum sensitivity to geometric distortion, minimum number of independent 
stress parameters, and rank sufficiency, (3) shear-deformed shell finite 
elements for laminated composites based on the total Lagrangian description 
accounting for anisotropic material behavior, dynamic response, arbitrary 
laminate configuration, and arbitrary ply properties, and (4) large-aspect- 
ratio finite elements for nonlinear shell-type structures analysis based on a 
higher order "degenerated: shell element with nine nodes and accounting for 
elastic-plastic behavior. 

Algorithms /convergence include (1) self-adaptive solution strategies 
focusing on alternative formulations for developing algorithms that avoid the 
need for global updating and inversion, (2) element-by-element solution proce- 
dures where approximate factorization is considered for solving the large 
finite-element equation systems that arise in nonlinear structural analyses, 

(3) automatic finite-element generators, where the use of VAXIMA (MACSIMA in 
VAX) is used to generate the equations that describe the finite-element formu- 
lation, and (4) convergence criteria considering the effects of underintegra- 
tion on the element approximations. 

Inelastic analysis and interactive elements include (1) inelastic and 
dynamic fracture, focusing on large, dynamic crack propagation by developing 
methods for new path-independent integrals, for generalized inelastic consti- 
tutive relations, and for new complementary energy approaches, (2) interactive 
finite elements where novel methods are developed to describe the bearing/- 
support interaction for general engine structural dynamics analysis, (3) non- 
linear composite structures analysis applicable to high-temperature composite 
material behavior, and (4) three-dimensional inelastic analysis via boundary 
integral, concentrating on the development of discrete element analysis based 
on boundary integral concepts having the potential for substantial computa- 
tional efficiency for local concentrations compared with traditional finite- 
element analysis. 


i i i 



Collectively, the papers included in these proceedings are representative 
of cutting-edge methodology in all three disciplines. The authors of the 
papers are nationally and internationally recognized experts in their respec- 
tive areas. The reader should bear in mind that the papers describe research 
in progress. Results and conclusions reported are subject to revision as 
additional results become available. In any event, these results and conclu- 
sions are those of the respective authors and not of the U.S. Government. 

As workshop coordinator, I would like to thank all the authors for pre- 
paring their papers on time as well as the attendees whose extensive 
contributions to the workshop discussions helped make the workshop a very 
successful technical information exchange forum. Finally, I thank all those 
who helped with the mechanics of organizing and conducting the workshop. 


C. C. Chamis 

Lewis Research Center 

Cleveland, Ohio 
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SLAVE FINITE ELEMENTS: THE TEMPORAL ELEMENT 


APPROACH TO NONLINEAR ANALYSIS* 

Slade Gellin 
Bell Aerospace Textron 


SUMMARY 


A formulation method for finite elements in space and time incorporating non- 
linear geometric and material behavior is presented. The method uses interpolation 
polynomials for approximating the behavior of various quantities over the element 
domain, and only explicit integration over space and time. While applications are 
general, the plate and shell elements that are currently being programmed will be 
used to model turbine blades, vanes and cumbustor liners. 


INTRODUCTION 


The extension of the finite element method to the solution of transient field 
problems by discretizing in time as well as in space has been investigated in the 
last fifteen years. The works most often cited are those by Argyris and Scharpf 
(Ref. 1) and Fried (Ref. 2). The starting point of their work was Hamilton's 
principle. Generally, sample problems consisted of axial thrust members, and only 
linear geometric and material properties were assumed. There were some arbitrary 
decisions made concerning the proper use of the initial conditions in the global 
system in order to derive the correct set of equations for the given problem. Many 
of the theoretical fine points were explained by Simkins (Ref. 3), who demonstrated 
his findings using point mass structures. His results indicated that very high 
levels of accuracy can be obtained with only a very small number of elements. He 
notes that the temporal elements are well suited to handle sudden changes in load 
function, extending the interval of solution indefinitely without restart, and 
providing great detail to the solution in any subinterval. Furthermore, conven- 
tional step-by— step integration algorithms may call for a large number of time 
steps, particularly for the hyperbolic equations of structural dynamics should the 
excitation or material properties change rapidly in time. It is within this spirit 
that this research was undertaken. 

During the course of this work it was found that many theoretical questions 
needed to be answered. Some of these questions have deep physical and mathematical 
importance, others are exercises in intellectual gamesmanship. Hopefully, they will 
be addressed in forth-coming papers. The focus here will be on incorporating non- 
linear geometric and material behavior into the temporal element approach, within 
the framework of the ground rules discussed in the following paragraphs. 


*This work was performed under NASA Contract No. NAS3-23279. 
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First, the elements will be of what will be referred to as the "prismatic" 
type. The space time domain is thought of as a long, prismatic, member in four 
dimensions, with the "length" dimension corresponding to time. The three dimen- 
sional "cross-section" is the spatial model of the structure. A prismatic element 
consists of the space-time domain occupied by a spatial finite element from t=t-^ 
to t=t j . Geometrically, this is a well defined set of "longitudinal fibers" over 
a certain "length". A set of local coordinates is chosen with the time direction 
the same as global time, but translated so that the time interval goes from 0 to 
T=tj-t i . See Figure 1. 

Second, all field quantities will either be derived from the displacement 
field directly or calculated at grid points and then interpolated with their own 
shape functions over the space-time domain. This will facilitate the explicit 
integration that also will be required for the procedure. Note that certain 
spatial geometries for the element may require special handling. Isoparametric 
representation of elements may not be desirable as only numerical integration 
schemes will generally be feasible. 

Finally, it will be assumed that the constitutive laws will be linear be- 
tween the time derivatives of the appropriate kinematical and dynamical quantities. 
Some linearization must be available in order for the matrix techniques of finite 
element analysis to be applicable. 

The non-linear algorithm is derived in the following section. In the sub- 
sequent section, a simple, but purposely poorly constructed example conducive to 
hand computation is presented. Finally, some brief comments about the current 
research being done using the temporal element concept will be given. 


Non-Linear Algorithm 


Like most non-linear algorithms, the one presented here is based on an 
iterative procedure involving quasi-linearization. The difference in philosophy 
between the method presented here and the conventional step-by-step methods, such 
as the tangent modulus or residual force methods, is demonstrated graphically in 
Figure 2. In both Figures 2(a) and 2(b), the heavily drawn curve represents the 
"exact" time history of quantity A, which may be thought of as the displacement 
of a certain point, or a stress at a point, etc. In Figure 2(a), the step-by-step 
method has calculated the path OA, representing the history up to that point. At 
A, the procedure generates successive approximations ABi until convergence to the 
path AB is achieved. The procedure is now repeated at B. For the algorithm 
presented here. Figure 2(b) indicates that successive iterations generate an 
entire load history. This load history allows for determining loading and unloading 
paths for the entire time interval of interest, eliminating the guesswork or sub- 
interval changes usually associated with the step-by-step methods. Iterations 
converge until the curve OC is obtained. 

The theoretical basis for the algorithm is Hamilton’s law of varying action. 
Essentially, it is the principle of virtual work integrated over time. It is 
expressed mathematically as 
T 

S f (<$£ a - 6v p) dVdt - {6u}^{F^} = 0 (1) 

°1 V 
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where e, a, v and p are the strain, stress, velocity and momentum fields, 
respectively, and {u} is the set of element degrees of freedom (d.o.f.). Note that 
each element of {u} represents a displacement measure at a certain point in place 
and time. The vector {F' ? } represents the equivalent point impulse-momentum 
difference due to known body forces and tractions applied over S a, the portion of 
the surface S where tractions are specified and thermal loading as well as "known 
momenta" applied over the time boundaries and equivalent loads exerted on the 
element by its contiguous neighbors. The variations in (1) are taken with respect 
to the displacement field "u, which is expressed in terms of shape functions 
[N(x,t)] as 

u = [N]{u} (2 ) 

The displacement u is assumed to be admissible. In this case, that means the shape 
functions [N] are interelement continuous, as well as continuously differentiable 
in the space-time domain. The strain-displacement law may be written formally as 

e = f(u) ( 3 ) 

■ ^ 

where f is a function of u and its spatial partial derivatives, and, in general, 
is non-linear. Taking variations of (3) yields 


8e = f''(u) 6 u ( 4 ) 

where the prime has a general meaning relating to derivatives with respect to vT and 
its partial derivatives. In a similar manner, the velocity— displacement relations 
may be expressed as 


v = g(u) ( 5 ) 

where g, like f, is a (non-linear) function of u and its spatial and temporal first 
derivatives. Taking the variation of (5) yields 


<Sv = g ^(u)<5u 

It is also interesting to note the time derivatives of ( 3 ) and (5) 


£ = f'(u) u 

o ^ q 

V = g * (u) u 


( 6 ) 


(7a) 


(7b) 


Constitutive laws between £ and a, and v and p may be formulated, according to the 
groundrules specified above, as 


o. ^ _ x 0 ^ 

a = [E(a,e,0,x,t, .. .)]eT+ t t,...) 

o ^ 

P = [p(a,e,0,x,t, . . ,)]v“+ if (a,£r,0,x,t, . . .) 


( 8 a) 

( 8 b) 
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where 0 is the temperature field. Equations (8) are integrated from 0 to t; thus, 


t ° 

a = f o dt + a (9a) 

o 0 

c 9. 

p = / p dt" + p (9b) 

o . ° 

where a and p are the values of the stress and momentum at local time t=0. These 
may be approximated by 

a Q = [N 0 (x)](a o } (10a) 

P Q = [N p (x)Ip o > (10b) 

where the number of parameters in {a 0 } and {p } is arbitrary. Ideally, [N ] is 
interelement continuous, and both [N^] and [N°] are as sophisticated as thl stress 
and momentum fields derived from [N] in the linear theory, though neither require- 
ment is really mandatory. 

In the iterative procedure to be used, "current" values of displacement, 
stress, strain, etc., for the entire load history are in hand. These values are 
used to evaluate ?" , g - ", [E], [p] and The quantities 5tT and are calculated 

from (2) where the set {u} are unknown, representing the "updated 11 solution. 
Similarly, {a Q } and {p Q } are unknowns. 

Equations (3) -(7) may be re-expressed as 


6e = 

[f "][N]{6u) 

(11a) 

o 

e = 

[f"][N]{u> 

(lib) 

6v = 

[g"][N]{6u} 

(11c) 

O 

v = 

Cg"][N](u} 

(lid) 


where [f"] and [g"] are 6x3 and 3x3 operator matrices, respectively. Stress and 
momentum matrices are defined by 
t o 

[S] = / [E][f"][N] dt" (12a) 

° t 

[M] = / [p][g"][N] dt" (12b) 


and stress and momentum vectors are defined as 
1 SL 

T = / T dt" 

_ °to ^ 
ir = / tt dt 
o 


(13a) 

(13b) 
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(14a) 


Using (10) , (12) and (13) in (9) yields 
a = [s]{u} + t + [n ]{a } 

Q O 

P= [M]{u} + v + [N p ]{p o ) (14b) 

Equations (14) may be evaluated at t=T. (Quantities evaluated at this time are 
given a T subscript.) Assuming that [N ] and [N ] can also approximate cf T and p" , 
equations (14) take on the form ^ T T 

[S T ] (u) + T t + [N a ] {u o -a T > = {0} ( 15 a) 

C m t ]{u) tt t + [N p ] {P 0 -P T > = {0} (15b) 

at t=T . Equations (15) are used as subsidiary conditions to the problem. They 
are added to (1) with the use of lagrange multiplier fields. In particular, the 
choices for the stress and momentum fields, respectively, is given as 


x a ■ IW 

(16a) 

* P ■ ^VV 

(16b) 


These fields are used with (15) and integrated over the volume. These conditions, 
as well as equations (11a), (11c), and (14) are used in (1) to give 

(6u } T ( [K ] {u } - {F*} - (F®}-{F^}+[A J]{a } - [A e ]{p }) 
c 1 Tl O O p o 

+ <$({X a ) T ([b®]{u) + [c^]{a Q -a T } - {q^} ) 

+U p } T ([B p ]{u} + [Cp]{p o -p T } - {q®}) = 0 ( 17 ) 

where 


T 

[K e ] = / / ([N] T [f"] T [s] - [N] T [g'] T [M] ) dVdt 

(18a) 

[A a ] = i i CN] T [f'] T [N a ]dVdt 

(18b) 

T 

Up] ' i !, [N] [s'f[N p ]dVdt 

(18c) 

[B a ] ■ l [»a ]T [ S T ]dv 

(18d) 
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[B®] = f [N p ] T [M T ]dV 


(18e) 


[cj] = i ^ N a ] [N a ]dV 


(18f) 


[C®] = / [N r[N ]dV 
P v P P 

T 

{F^} = -/ /[N] T [f^] T T dVdt 

T O V 


{F®} = / JtN] [g'] TTdVdt 

F O V 


{ qp=-/[N a ] X T dV 


{ q® } =-/[N p ] V V 


The set of equations generated by (17) when variations are taken on (u), {<?_,}> 

{p 0 }, and [^O are 


T T 

.G-I r «, *1 , r-n^l I 


[K e ](u} + [A^]{a o J - [Ap{p o )+[B“] <A a } + [Bp] U p J 


- {F*}+{F®}+{F*} ' 

- { 0 } 

[c e ] {X ) - {0} 1 

p p 

[B*](u] + [C^](a o -0 T } = {q^} 

[ B p](u}+ [Cp]{p 0 -P T )= (qj) 

Note that the matrices [C e ] and [Cp^will be square and invertible, thus making 
the mulipliers identically zero. They may then be omitted from equation (19a). 
Equations (19a, d, e) are assembled into the systems 

[K] {u } + [A a ](o} - [Ap] (p } = (?) 

[B a ]{u] + [C a ] (a) = {q T > 


[B p ] (u ) + [ c p Hp 0 ) = 



where the assembly enforces the conditions of continuity across a time boundary 
for stress and continuity of momentum across a time boundary to the extent that no 
impulses per unit volume are applied, and, if such impulses are implied, an 
appropriate discontinuity is maintained. The results should leave [C a ] and [C ] 
square and invertible. Thus, equations (20b, c) are solved for {a} and {p } and P 
used in (20a) to derive the global stiffness equations 

[K] (u } = {p} (21) 

where 

W - M - [A 0 ] [C a ]' 1 tB a ] + [A p ][C p ] _1 [B p ] (22a) 

[P] - [P] - [A 0 ] [C 0 ]' 1 [q T ] + [AplCCp]" 1 ^] (22b) 

Equations (21) are solved for {u} and then back-substituted into all the pertinent 
equations to calculate the important quantities to be used in the next iteration. 


A Simple Numerical Example 


To demonstrate some of the procedures developed, a numerical example will 
be worked out, A rod of length L, and cross-sectional area A built in at both 
ends, is loaded at its midpoint by a continuously time varying load p(t) given as 


?(t) = 2a A 
Y 



(23) 


where T is the time interval of interest and is the nominal yield stress of the 
material satisfying the uniaxial stress strain law of the Ramberg-Osgood type 
(Reference 4) 



where E is the elastic modulus of the material. The loading is assumed quasistatic 
and thus dynamic effects may be ignored. (Note that if P(t) is discontinuous, 
dynamic effects must be introduced in order to maintain continuity of displacement 
and stress across a time boundary.) The rod is initially in the undeformed state; 
thus, not only is u(x,0)=0, but a(x,0)=0. See Figure 3. 

The problem is discretized using four elements, each of dimensions L/2 by 
T/2, as seen in Figure 4. Each element is of the type shown in Figure 5. The 
displacement field for this element is modeled using the bilinear shape functions 

u(x,t)=u® 1 (l-x/L')(l-t/T , )+U 2 1 d-x/L'Mt/T') 

+ u® 2 (x/L 1 ) (1-t/T* )+ u 22 (x/L ' ) (t/T ' ) (25) 
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where the e superscript is a reminder that the quantities in (25) refer to local 
numbering and local axes. It can be shown, using methods derived in (Ref. 3), 
that the element statics equations for a linear elastic rod are 



(26) 


For the structure being studied here, L'= L/2 and T'= T/2; furthermore, only U £2 
and U 32 are non-zero. The global stiffness equations for the linear elastic 
system are thus 



(27) 


The solution is U 22 = 1/4 a y L/E, u 32 = 1/2 a y L/E. The stress state is a 2 = c(T/2) 
= l/2a Y and 03 = a(T) = a Y , where the values are in tension on the left half of 
the rod and in compression for those on the right. The results are exact. 


When the non-linear behavior is incorporated, the formulation described here- 
in is used. The non-trivial matrices to be found are [K e ], [ J , [B 0 ] and [C®]. 

The rate dependent form of (24) is 

a = E a e (28) 


where 

a(t) 



(29) 


For the problem at hand, the stress state is assumed (for each element) to be a 
function of time only. Thus, 

[N 0 ] = 1 (30) 

The matrix [E] is approximated as a linear function in time by interpolating through 
the end points of the time interval; specifically, 

[E] = Ea(0) [1 - t/T'] + Ea(0(t/T") (31) 

Note how extremely poor this approximation is for the desired time step used 
in this example. 

For the bilinear shape functions, 

[f'][ N] = [1 -1 -1 1] (32) 


8 



Incorporating (31) and (32) into (12a) yields 


[«- [(f)- 1 (§') 2 ] 


q(T ) / t 


[i -l -l l] 


Evaluating (33) at t=T , noting (30), and integrating over the volume yields 
[B®] = EA [! _! _i i] 


[C®] = AL' 


Using (33) in (18a) gives 


[K e ] = 


9ki. + qq 

\ 8 2 V 

feq a + q_A 

\ 24 a) 

So. + ^L 

8 24 

5a 0 ai_ 

r\ / * n 


Sq. . aj_ 

n * * a / 


5a Q ai 

24 8 


ia + 

8 24 


5a 0 ai 

24 8 


qp ■ qi 
8 24 


5q Q , oti 

24 8 


q Q + ai 


5a 0 ai 

24 8 


ao ai 

A * A # 


5a 0 ai_ 

24 8, 


qp , qi 
8 24 


5a 0 , ai 

24 8 . 


where a Q = a(0) and ai 


= a(T^) , and using (30) in (18b) yields 


[A®] = AT' 


Global assembly is now performed. First, it is remembered that only u „2 and u„„ 
are non-trivial; secondly, it is known that cr(0) = Oi = 0 and thus a(0j = ai = I. 
Let a 2 subscript be used for quantities associated with t=T/2, and a 3 subscript 
be associated with quantities evaluated at t=T. Furhter, by symmetry, the value 
of the tension in an element on the left hand side of the structure would equal the 
value of compression in the corresponding element on the righ-hand side of the 
structure, so = ^Right at a gi ven time t. Finally, it is to be remembered 

that for the current discretization T'’ = T/2 and L' = L/2. All these imply that 
(20a) and (20b) may be written as 


d3 

q 2 

+ 

q3 " 

j u 22 

12 

4 

12 

q 3 

5a 2 

+ 

a3 

1 

4' 

12 

4"J 

\ 32 J 


+ AT 


0-2 = 0 AT 


j 1/2 ( 

( 5 /I 2 ! 
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mh+ -2 z> [10] | “22 1 -f 02 - 0 (37b) 

Solving (37b) for C 2 and using the results in (37a) yields the global stiffness 
equations 



f 11 

+ £2 

as 

a 2 a 3 

( u 22 J 

i ( 1/2 i 

| 


12 

2 

" 12 

4 12 

i 

I 

AET 






It 

(38) 

L 

1 

+ ^- 

a 3 

5a 2 , a 3 

r 32 ! 

1 5/12 | 

' 


_ 2 

12 

4 

12 4 . 

1 \ j 



Note that if the structure is considered linear elastic that a 2 = CX 3 = 1, and 
equations (38) become identical to the elastic global equations (27). 

Equations (38) were solved iteratively for u 22 9 ° 2 anc * ° 3 us ^- n § a 

hand-held calculator maintaining four -place accuracy in the coefficients of the 
stiffness matrix. Table 1 lists the results of each iteration where the elastic 
case is taken as the initial condition. It should be noted that the trend indi- 
cated by the table implies that greater accuracy must be retained to achieve con- 
vergence; however, three place accuracy is reached rather quickly. The results 
are not bad when considering the absolutely horrendous approximations used in the 
example. It is concluded that the general method is viable. 


Related Research 


As noted in the introduction, there were many theoretical questions which 
were needed to be answered that arose during the course of this study. Many of the 
questions concerning the incorporation of the initial conditions to the transient 
problem into a formulation that is inherently suited to boundary type conditions 
were answered both independently and through (ref. 3) . There were other questions 
concerning non-prismatic or semi -prismatic discretizations, and their possible use 
in both formulating elements and in modeling specific problems in space and time. 
Several modified or hybrid variational formulations were developed to model ele- 
ments possessing certain qualities, such as an embedded hole or non-rectangular 
geometry for grid elements, where interelement continuity would be more difficult 
to maintain. A search for a four -dimens ional , unified theory for elastic dynamic 
solid mechanics was undertaken in order to better understand the parallel conditions 
that must exist between the spatial and temporal properties of a structure. Not to 
be overlooked are the various computer programming difficulties, particularly the 
large numbers of d.o.f. per element that could be generated for elements of this 
type, and the adaptability of element modules with a solution procedure provided 
by another host program. 


CONCLUDING REMARKS 


The Slave Finite Element approach to non-linear analysis has been presented. 
Despite the large number of dof in the system for any given problem, the anticipated 
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reduced number of calculations, the ease in changing time steps by using various 
grid options, the increased detail within any subinterval, and the accuracy demon- 
strated in elastic problems demonstrates the viability of the method. Continued 
research is necessary to fully understand and exploit the potential of the method. 
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TABLE 1: NUMERICAL COMPUTATIONS FOR SAMPLE PROBLEM 


Iteration # 

U2 aE/ a L 

U3 2E/a„L 

a 2 /c y 

< 73 /a y 

02 

(*3 

0 

.2500 

.5000 

.5000 

1.0000 

.9852 

.2059 

1 

.2409 

.6252 

.4782 

.9359 

.9896 

.3058 

2 

.2423 

.6016 

.4821 

.9475 

.9889 

.2853 

3 

.2420 

.6061 

.4813 

.9452 

.9890 

.2892 

4 

.2420 

.6053 

.4813 

.9457 

.9890 

.2884 

5 

.2420 

.6054 

.4813 

.9457 

.9890 

.2887 



(no further 

change in stiffness 

matrix) 



EXACT 

.2504 

.7143 

.5000 

1.0000 



% ERROR 

3.4 

15.3 

3.7 

5.5 
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(a) 



Figure 2. Comparison of Non-Linear Algorithms (a) Step-By-Step Method 
(b) Slave Finite Element Method 
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NEW VARIATIONAL FORMULATIONS OF HYBRID STRESS ELEMENTS* 

T.H.H. Pian,^ K. Sumihara^ and D. Kang”^ 
Massachusetts Institute of Technology 


SUMMARY 


In the variational formulations of finite elements by the Hu-Washizu and 
He! 1 inger-Reissner principles the stress equilibrium condition can be maintained 
by the inclusion of internal displacements which function as the Lagrange multi- 
pliers for the constraints. These new versions permit the use of natural coordi- 
nates and the relaxation of the equilibrium conditions and, render considerable 
improvements in the assumed stress hybrid elements. These include the derivation 
of invariant hybrid elements which possess the ideal qualities such as minimum 
sensitivity to geometric distortions, minimum number of independent stress param- 
eters, rank sufficient and ability to represent constant strain-states and bending 
moments. Another application is the formulation of semiLoof thin shell elements 
which can yield excellent results for many severe test cases because the rigid 
body nodes, the momentless membrane strains and the inextensional bending modes 
can all be represented. 


INTRODUCTION 


In the original formulation of the assumed stress hybrid elements by the 
modified complementary energy principle (ref. 1) and the later extension using 
the Hel 1 inger-Reissner (ref. 2) the assumed stress are made to satisfy the equilib- 
rium equations a priori. This restricts the derivation to the use of physical 
coordinates such as Cartesian coordinates and shell surface coordinates, etc. and 
to the coupling of the various stress components by the equilibrium equations. 

Such restrictions are the chief obstacles for the construction of shell elements 
by the hybrid method and the main reasons why the element properties often deterio- 
rate rapidly when unfavorable reference coordinates are used for stresses and/or 
when the element geometry is distorted. 

Recently new versions of the Hell inger-Reissner and the Hu-Washizu principles 
have been suggested for more efficient formulations of hybrid stress elements 
(refs. 3 and 4). This paper presents two specific applications of these principles 
(1) the formulation of invariant hybrid stress elements by the Hell inger-Reissner 
principle and (2) the derivation of semiLoof shell elements by the Hu-Washizu 
principle. 

The research reported here is one part of a program the objective of which 
is to advance the numerical tools for analyzing engine structures which needs to 
be modelled by 3-D solids and/or shell structures. 


★ 

The work was sponsored by the NASA Lewis Research Center under NASA Grant 
No. NAG 3-33. 

4 * 

Professor of Aeronautics and Astronautics 
^Graduate Student 
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FORMULATION OF INVARIANT HYBRID STRESS ELEMENTS 


The modified version of the Hellinger-Reissner principle for the formulation 
of the element stiffness matrix may be stated as 


TTr 


= / 
V 


r 1 rJc 
lr 2 2 ~ 


a + cJ(D u„) - (D^a)^ u, ]dV = Stationary 
~ - ~ ~q ~ ~ ~a 


( 1 ) 


where a = stress, S = elastic compliance, V = element volume, Ug = element 
displacements that are interpolated compatibly in terms of nodal displacements q 
and u^ = additional internal displacements which are expressed in terms of parame- 
ters" X that can be statically condensed. Here 

(D T a) = 0 (2) 


is the equation of equilibrium. Hence the last term in the functional plays the 
role of conditions of constraint with the corresponding Lagrange multipliers. Thus, 
the assumed stresses g need not be coupled initially and the introduction of dis- 
placement parameters X will reduce the number of independent stress parameters. 


As an example, the derivation of the quadrilateral membrane element is 
presented. Here, for the a stresses, the nine terms are complete in linear terms 
in the natural or isoparametric coordinates £, n. 


1 

a = i 

f a 
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* — 
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The element displacements are interpolated by bilinear polynomials in the 
natural coordinates. Here it is expected that in the limiting case of rectangular 
elements the present formulation should yield the 5-B hybrid stress element (ref. 5) 
that is known to have excellent properties. 

Thus, four internal displacement terms are need. They are, indeed, the same 
ones used by Wilson et al . (ref. 6) in their incompatible elements, i.e. 


u x = X-,0-5 2 ) + X 2 (l-n 2 ) 

v x = x 3 (l-c 2 ) + x 4 (l-n 2 ) 


(4) 


It is noted that these additional displacement terms are essential because they now 
make the displacements u complete in quadratic terms. 

In the actual formulation it turns out that the four X terms yield only two 
independent equations of constraints for the B's. It was found necessary to 
consider the element geometry with a small perturbation as shown in Figure 1. With 
the x and y components of the perturbation represented by ±Ax and ±Ay the Jacobian 
is 
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( 5 ) 


where 



a, + 

a 2 n 

+ Ax (1 ■ 

-n 2 ) 


_ a 3 + 

a 2? 

-2Ax£n 


_ 1 
4 

(“X-j ^ 

h x 2 

+ X 3 " 

X4) ; 

_ 1 
4 

■ < X 1 - 

x 2 

+ x 3 " 

x 4 ) ’> 

1 

4 

■ ( - x l 

- x 2 

+ x 3 + 

x 4) > 


b-| + b 2 n + Ay(l-n 2 ) 

b 3 + b 2^ " 2 ^ n 

b l = + *2 + *3 ' y 4> 

b 2 = 1 <*i - *2 * H ' H > 

b 3 = ? ("J'l ' *2 + *3 + V 


( 6 ) 


Here x.,y. are the coordinates of the corner nodes and the ratio Ay/Ax is equal 
to b-j/a-j . 


In this case the four X-terms yield four equations of constraints and 8-, 8^» 
8 8 and Bo can be expressed in terms of 8 q and e_. The assumed stresses can be 6 
rearranged in matrix form as o 5 


f 1 

a x 


0 

1 

0 

a 2 n 

a 3 ? 


1 

?1 

a 

y 

r ~ 

0 

1 

0 

b 2 n 

b 3^ 

' 

* 

T 

( x y J 

1 

0 

0 

1 

a-j^n 

a 3 b 3 C _ 

1 

1 

Be 
5 ; 


(7) 


In the case of a rectangular element with £ and n in parallel with x and y, 
b-j and a^ vanish and the 5 8-terms of reference 5 is obtained. For hybrid stress 
elements, particularly of higher orders, there exist many possible choices for the 
assumed stress. The new variational principle, thus, provides a rational procedure 
to establish the appropriate stress terms for hybrid elements both of regular 
shape and with geometrical distortions. Because the present formulation is based 
on natural coordinates the resultant element stiffness matrix is always an 
invariant. It can also be shown that it always has sufficient rank. The element, 
of course, can also pass the constant strain patch test. 


It turns out that the same result can be obtained when the stresses based 
upon the basis vectors of the natural coordinates £ and n are expressed as 


T = \ 


11 1 
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0 
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(8) 


and in converting the tensor stress x 1J to the physical component a 1J by 


a lj = J l J l r k£ 


(9) 


the value of the Jacobian at (^,n) = (0,0) is used. By such a step the constant 
stress state can be maintained hence the resulting element can pass the patch test. 
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A tapered and swept panel with one edge clamped and the opposite edge acted by 
a distributed shear load was used bt Cook (ref. 7) for testing the sensitivities of 
finite elements to geometric distortions. The resulting normal stress distributions 
along the mid-span section are determined by the following four elements using the 
same 4x4 mesh as shown in figure 2: 

(1) Element Q-4 by bi -linear assumed displacements. 

(2) Element HG by original 5-3 assumed stress method with the global x-y 
axes as reference. 

(3) Element HL by original 5-3 assumed stress method with the reference 
x'-y' axes located at equal angles with the natural £-n axes. 

(4) Present invariant hybrid stress element. 

A solution obtained by element HL using 16x16 mesh is used as a reference. 

It is seen by the result by the present element is already very close to the 
reference solution while the solutions by the other elements all have much 
larger errors. The assumed stress given by eq. (8) is thus the optimal choice 
for assumed stress hybrid elements. 

Eight-node solid elements have also been constructed by the same approach. It 
is used to analyze the bending of a rectangular bar with two elements which are 
distorted. The effect of the element distortion on the tip deflection and normal 
stress is shown in figures 3 and 4. Also included for comparison are the results 
obtained by the assumed displacement method and by the original hybrid stress 
method using the beam axis as reference. It is seen that under large distortions 
the original hybrid stress element may become almost as rigid as the assumed 
displacement element, but the present invariant element is much less sensitive 
to distortions. 


THIN SHELL ANALYSES BY SEMILOOF ELEMENTS 

Thin shell elements are formulated by the Hu-Washizu principle with additional 
internal displacements for enforcing the stress equilibrium conditions and with 
relaxed continuity conditions along the interelement boundaries. The principle 
is expressed as 


tt hw = / [ | e T c e - o T e + o T (D u q ) - (D T a) T u x ]dV - / T T (u q -u)dS 

= stationary (10) 

where strains e, stresses a and element displacements Uq and u^ and boundary 
displacement ITare independent. The key step in the finite element formulation 
is that both ~a and e are approximated by the same function, i.e. 

e = Pa and a = P3 (11 ) 

where P are not coupled among different components. The displacements u , u^ and 
u are "interpolated by 
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( 12 ) 


u = Nq ; u, = MX ; and u = Lq 

~M — ~A 

and the functional tt^ is reduced to 


where 


^HW = \ a " ^ a + ^ 

J = / P T C P dV ; H = / 

~ ' V ~ V 

G = / P T (D N)dV ; R = / 

~ V ~ ~ ~ ~ V 


q - 6 T R X 
P T P dV 

(D T P) T M dV 


(13) 


(14) 


Here H is a symmetric and positive definite matrix which can be partitioned with 
submatrics located only along the diagonal. Thus, the inversion of H reduces to 
that of the individual submatrices. 


Variation of tthu with respect to 3, a and X in the element level will enable 
the solutions of these variables in terms of q and the following expression can be 
obtained for the element stiffness. 

k = G T W G - G T W R (R T W R) _1 R T W G (15) 


where 


W = H _1 J H” 1 (16) 

and the order of R* W R is the same as the number of equilibrium constraint equations 
obtained by X. 

To maintain good qualities a thin shell element must be able to represent 
(1) the rigid body modes, (2) the momentless membrane strains and (3) the inexten- 
sional bending modes. By using the Hu-Washizu principle for' which the assumed 
strain and stress components are uncoupled the last two deformation modes can be 
easily included and by satisfying the equilibrium conditions in the element level 
the rigid body displacements can be guaranteed. 

Another condition that is desirable for the formulation of thin shell elements 
is for the three displacement components u, v, and w to be of the same order. For 
ordinary shell elements for which the nodal displacements w and w. and w . are 
needed to interpolate the lateral displacements while only u and ’ a v are ’^used for 
the assumed inplane displacements, these displacement components will naturally be 
unbalanced. In a semiLoof element, however, all the displacement components are 
interpolated in terms of their own nodal values hence are automatically of identical 
orders. By the assumed stress hybrid method the continuation of the normal rota- 
tions w >n along the boundary are imposed through the boundary displacements u. 

Thus, in addition to the advantage of expressing equilibrium equations at the 
boundary nodes with only three degrees of freedom, a semiLoof element is also most 
natural for the balanced assumed displacements. 

Hybrid semiLoof elements for triangular and quadrilateral planform with 


21 



24 and 32 DOF respectively were constructed based on the shallow shell theory of 
Marguerre and the optimal number of p's and A's were determined. The nodal displace- 
ments are u, v and w at each corner node and each mid-side node and w at 1/3 way 
points of all edges. ,n 

The semiLoof approach has also been used to construct two hybrid quadrilateral 
elements of 32- DOF based on the cylindrical shell theory using 36 and 38 6-parame- 
ters respectively. The stresses are later constrained by using 7A's. For the 
36-6 element the in-plane and moment components are complete in quadratic terms. 

The element has two kinematic deformation modes, hence a 38-6 element is also 
developed by adding cubic terms for the in-plane stresses components 5,^ and 5,2 2 # 

For the internal displacements the u^ and u« components each contains 3 terms 
while the w component contains only a single term. 

Figures 5. and 6 .present comparisons of vertical deflections at the edge of a 
cylindrical shell roof under gravity loads by three 32 DOF semiLoof quadrilateral 
elements, one which is based on the shallow shell theory and the other two, by the 
deep shell theory. Also included for comparison are the solutions by the assumed 
displacement methods. The result by Dawe's (ref. 8) is by.54-D0F triangular elements 
based on the deep shallow theory. Dawe's element is known to be an excellent element 
although it is not easily adopted because of its use of higher derivatives as nodal 
displacements and its difficulty for handling the intersection of shells. The 
result by Cowper et al . (ref. 9) is again obtained by the triangular elements based 
on the shallow shell theory. That element as well as the degenerated solid element 
with reduced integration (ref. 10) all coverge to that by the shallow shell theory. 

On the other hand, the solutions by the hybrid semiLoof elements based on the 
cylindrical shell theory, converge to the solution obtained by the deep shell 
theory. This is a clear indication that the appropriate shell theory must be 
incorporated in the construction of shell elements. Figure 7 shows stress distri- 
butions obtained by the semiLoof element with 386 using 3x3 and 4x4 meshes. 

Morley (ref. 11) has proposed a challenging test example which consists of 
a thin circular cylinder with reciprocal shear tractions applied along a cut in the 
axial direction resulting in a state of uniform torsion (fig. 8). This problem is 
very sensitive to the assumed displacement fields in the finite element. Another 
severe test case is a thin cylindrical shell under internal pressure for which the 
element used must be able to represent momentless membrane stresses. It has 
been found that the present semiLoof cylindrical shell element with 366 and 7A 
can pass these tests successfully. 


CONCLUSIONS 


The new versions of the Hellinger-Reissner and Hu-Washizu principles using 
additional internal displacements open up the options for constructing hybrid 
stress elements. A new hybrid stress element developed by the present method has 
invariant properties and has been demonstrated to be much less sensitive to 
geometric distortions. Several semiLoof shell elements have been constructed and 
proved to be successful even for very severe test examples. Introduction of the 
new variational principles, thus, represents a significant advancement in the 
finite element development. 
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A SHEAR DEFORMABLE SHELL ELEMENT 
FOR LAMINATED COMPOSITES^ 


W. C. Chao and J. N. Reddy 

Department of Engineering Science and Mechanics 
Virginia Polytechnic Institute and State University 
Blacksburg, VA 24061 

SUMMARY 

A three-dimensional element based on the total Lagrangian description of the 
motion of a layered anisotropic composite mediim is developed, validated, and used 
to analyze layered composite shells. The element contains the following features: 
geometric nonlinearity, dynamic (transient) behavior, and arbitrary lamination 
scheme and lamina properties. Numerical results of nonlinear bending, natural 
vibration, and transient response are presented to illustrate the capabilities of 
the element. 


INTRODUCTION 

Composite materials and reinforced plastics are increasingly used in 
automobiles, aircrafts, space vehicles, and pressure vessels. With the increased use 
of fiber-reinforced composites as structural elements, studies involving the 
thermomechanical behavior of shell components made of composites are receiving 
considerable attention. Functional requirements and economic considerations of 
design have forced designers to use accurate but economical methods of determining 
stresses, natural frequencies, buckling loads, etc. The majority of the research 
papers in the open literature on shells is concerned with bending, vibration, and 
buckling of isotropic shells. As composite materials are making their way into many 
engineering structures, analyses of shells made of such materials become 
important. The application of advanced fiber composites in jet engine fan or 
compressor blades and high performance aircraft require studies involving transient 
response of composite shell structures to assess the capability of these materials 
under dynamic loads. 

A review of the literature indicates that first, there does not exist any 
finite-element analysis of geometrically nonlinear transient response of laminated 
anisotropic shells, and second, the 3-D degenerated element is not exploited for 
geometrically nonlinear analysis of laminated anisotropic shells. The present study 
was undertaken to develop a finite-element analysis capability for the static and 
dynamic analysis of geometrically nonlinear theory of layered anisotropic shells. A 


A more detailed account of this paper can be found in NASA CR-168182. 
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3-D degenerated element with total Lagrangian description is developed and used to 
analyze various shell problems. 

INCREMENTAL, TOTAL-LAG RANG IAN FORMULATION OF A CONTINUOUS MEDIUM 

Hie primary objective of this section is to review the formulation of equations 
governing geometrically nonlinear motion of a continuous median. In the interest of 
brevity only necessary equations are presented (see [1-5]). 


We describe the motion of a continuous body in a Cartesian coordinate system. 
The simultaneous position of all material points (i.e. , the configuration) of the 

body at time t is denoted by C t , and C q and C t+ ^ t denote the configurations at 

reference time t = 0 and time t + 6t, respectively. In the total Lagrangian 
description all dependent variables are referred to the reference configuration. 

t t t t 

The coordinates of a typical point in C t is denoted by x = ( x^ , x^ y x^). The 
displacement of a particle at time t is given by 


t t o t to 

u = x - x or u. = x. - x, 

~ ~ ~ 1 i i 


The increment of displacement during time t to t + fit is defined by 


u i = 


t4fit 


u i “ 


(i) 


( 2 ) 


The principle of virtual displacements can be employed to write the equilibrion 
equations at any fixed time t. The principle, applied to the large-displacements 
case, can be expressed mathematically as 


J 


V 


t+6t" . JTr , r t+5t„ 

u. 6u. dV + J S 

1 1 ° \ 
o 


. . 6( t+6t E . )dV 
ij ij o 


= J 


t+6t 


T. 6u, dA + / 
i i o J T 


t+fit 


F. 6u. dV 
1 1 o 


(3) 


o o 

where summation on repeated indices is implied; V Q , A Q , and p Q denote, respectively, 

a volume element, area element, and density in the initial configuration, are 

the components of second the PLola-KLrchhof f stress tensor, the components of 

the Green-Lagrangian stran tensor, T^ the components of boundary stresses, and F^ 

are the components of the body force vector. The superposed dots on u^ denote 

differentiation with respect to time, and 6 denotes the variational symbol. In 

writing Eq . (3) it is assuned that e.. is related to the displacement components by 
the kinematic relations ^ 


t+fit 1 , t+fit , t+6 t ^ t+6t 

ij 2 i,j j,i ra,i 


t+fit 


U .) 


where j = 5u^/ . 


(4) 
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The stress components 


can be decomposed into two parts: 

t+5t s = c s + s 

s lj s ij + s lj 


where is the incremental stress tensor. The incremental stress components s u 

are related to the incremental Gteen-Iagrange strain components, e.. = e. . + r) . . , by 
the generalized Hooke’s law: J J J 


S ij C ijkJl e kr 


ij ij 'll 
(6) 


where are the components of the elasticity tensor. Using Eqs. (4)-(6), one 

can be express Eq . (3) in the alternate form 

I p o t+5t u i 6u i dV o + I + I 'kl 6e i J )dV o 

o o 

+ / t s. . be.. dv = 6W - f t s.. 6 ti. . dV (7) 

\ ij ij o J v ij ij o 

v o o 

where SW is the virtual work due to external loads. 

FINITE-ELEMENT MODEL 

The coordinates of a typical point in the element can be written as (see 
Fig. 1) 


•¥ + 7, VW ¥ ‘^bottom < 8 > 


x i - , E , V 5 i’V-r (l ) t., + 7, 

J=1 J=1 


where n is the number of nodes, are the finite-element interpolation (or 

shape) functions, which, in the element take the value of unity at node i and zero 
at all other nodes, and ^ are the normalized curvilinear coordinates in the 

middle plane of the shell, and C is a linear coordinate in the thickness direction 

and xj, x^, and x^ are the global coordinates at node i. 

In the present study the current coordinates t x^ are interpolated by the 
expression 

C X = Z + | Ch -(9) 

1 j J 1 ^ J 

and the displacements by 


♦j ! M + { ( 

V u W Ch j 


°e J )] 
31/ ' 

(10) 

C e j )] 
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all four edges are clamped: 

u-v-w-* x -*«0 

6 . 8 . 10 . 12 . 

Center deflection, w (in mm) 

Load -deflect ion curve for a clamped 
cylindrical shell under uniform load 



Here t u^ and denote, respectively, the displacement and incremental displacement 
components in the x^-direction at the j-th node. 

Substitution of Eqs. (4)-(6) and ( 9)— (11) into Eq. (7) yields 


/ p o [T] t {u}dV o + ( C [K l ] + t [K NL ]){A} = t+6t {R} - t+6t { F} (12) 

V o 

where [K^] , { R} , and { f} are the linear and nonlinear stiffness matrices, 

force vector, and unbalanced force vectors 


Application of the Newnark direct integration scheme (see [2]) for the 
approximation of the time derivatives in Eq . (12) leads to 


[K] {A} = t+6t {R} - t46t { F } (k_1) 


(13) 


where 


[K] = a Q fc [M] + [K] 


t+6t 


{R} = t+St { R} + a 2 [ t {P 1 } -- ( t {P 2 } - t {P 3 })] + a 3 {P 4 } (14) 


and a Q , a£, etc. are the parameters in the Newmark integration scheme. 

DISCUSSION OF THE NIMERICAL RESULTS 

The results to be discussed are grouped into three major categories: 

(1) static bending, (2) natural vibration, and (3) transient response. All results, 
except for the vibrations, are presented in a graphical form. All of the results 
presented here were obtained on an IBM 37 0/3081 computer with double precision 
arithmatic. 

Static Analysis 

1. Cylindrical Shell Subjected to Radial Pressure Consider a circular 
cylindrical panel, clamped along all four edges and subjected to uniform radial 
inward pressure. The geometric and material properties are 

R = 2540 mm, a = b = 254 mm, h = 3.175 mm, 

9 = 0.1 rad, E = 3.10275 kK/mm 2 , v = 0.3 

Due to the symmetry of the geometry and deformation, only one quarter of the panel 

is analyzed. A load step of 0.5 KN/ra^ was used in order to get a close 
representation of the deformation path. Figure 2 shows the central deflection 
versus the pressure for the panel dimensions a = 254 mm and b = 254 mm. Hie 
solution for agrees very closely with that obtained by Dhatt [6] and the shell 
element of Reddy [7]. 
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2. Nine-Layer Cross-Ply (0°/9Q o /Q°/>>») Spherical Shell Subjected to Uniform 

Loading Consider a spherical shell laminated of nine layers of graphite-epoxy 
material = 40, Gj 2/ E 2 = G 13 = G 12 = G 23 f V 12 = subjected to 

uniformly distributed loading, and simply supported on all its edges (i.e. , 
transverse deflection and tangential rotations are zero). A comparison of the load- 
deflection curves obtained by the present element with those obtained by Noor [8] is 
presented (for the parameters h/a = 0.01 and R/a = 10) in Fig. 3. Ihe results agree 
very well with each other, the present 2-D results being closer to Noor's 
solution. This is expected because Noor f s element is based on a shell theory. 

3. Two-Layer Cross-Ply and Angle-Ply (45°/-45°) Shells Under Uniform Loading 

The geometry of the cylindrical shell used here is the same as that considered in 
Problem 1. The shell is assumed to be simply supported on all edges. The material 
properties of individual lamina are the same as those used in Problem 2. A mesh of 
2x2 nine-node elements in a quarter shell is used to model the problem. The results 
of the analysis are presented in the form of load-deflection curves in Fig. 4. From 
the results, one can conclude that the angle-ply shell is more stiffer than the 
cross-ply shell. The geometry and boundary conditions used for the spherical shells 
are the same as those used in Problem 2. The geometric parameters used are: R/a = 

10, a/h = 100. The load-deflection curves for the cross-ply and angle-ply shells 
are shown in Fig. 5. From the plots it is apparent that, for the load range 
considered, the angle-ply shell, being stiffer, does not exhibit much geometric 
nonlinearity. The load-deflection curve of the cross-ply shell exhibits varying 
degree of nonlinearity with the load. For load values between 100 and 150, the 
shell becomes relatively more flexible. 

Natural Vibration of Twisted Plates 


4. Natural vibration of cantilevered twisted plates Here we discuss the 
results obtained for natural frequencies of various twisted plates. This analysis 
was motivated by their relevance to natural vibrations of turbine blades. Consider 
an isotropic cylindrical panel with a twist angle 0 at the free end. Thble 1 
contains the natural frequencies of a square plate for various values of the twist 
angle 0 and ratios of side to thickness. A 2x2 mesh and 4x4 mesh of 9-node elements 
are employed to study the convergence trend. The results of the refined mesh are 
included in the parentheses. The results agree with many others published in a 
recent NASA report. 

Transient Analysis 


5. Spherical Cap Under Axisymmetric Pressure Loading Consider a spherical 
cap, clamped on the boundary and subjected to axisymmetric pressure loading, p Q . 

The geometric and material properties are 

R = 22.27 in., h = 0.41 in., E = 10.5 x 10^ psi, v = 0.3, 

p = 0.095 lb/in 3 , 9 = 26.67°, p Q = 100 ksi, 6t = 10 -5 sec. 

This problem has been analyzed by Stricklin, et al. [9] using an axisymmetric shell 
element. In the present study the spherical cap is discretized into five nine-node 
2-D or 3-D elements. Figure 6 shows the deflection of the center as a function of 
time. The present solutions are in excellent agreement in most places with that of 
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Figure 5. Center deflection versus load para- 
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Table 


Natural Frequencies of Twi sted Square Plates 

3 

(jj = coa 2 / ph/D , 0 = ^ i) ■ > v - 0* ** 3 

1 2(l-v ) 


a 

h 

Twist 

Angle 

1 

2 

Mode 

3 4 

5 

6 


0° 

* 3.4556 

+ (3.4583) 

8.4110 

(8.3353) 

22.0999 

(21.0238) 

28.2089 

(26.7465) 

31 .9740 
(30.1454) 

55.1625 

(52.0784) 


15° 

3.4359 

10.2920 

21.5199 

27.2054 

32.7430 

44.5375 

20 

30° 

3.3790 

(3.3694) 

13.7014 

(14.2222) 

19.9840 

(18.9795) 

25.0943 

(26.8104) 

34.3341 

(34.4591) 

45.8987 

(45.7547) 


45° 

3.2908 

18.1009 

15.9097 

23.5680 

35.5332 

45.7013 


60° 

6.1800 

17.8319 

15.5635 

24.1842 

36.1466 

44.9152 


0° 

* 3.33916 

**(3.3390) 

7.3948 

(7.3559) 

10.8083 

(10.883) 

18.4930 

(17.757) 

23.7907 

(22.769) 

26.0552 

(24.125) 


15° 

3.31713 

(3.3170) 

7.4816 

(7.4504 

10.8053 

(10.774) 

18.4043 

(17.771) 

23.6767 

(22.694) 

24.9474 

(24.083) 

5 

30° 

3.2538 

(3.2538) 

7.7593 

(7.7089) 

10.5248 

(10.478) 

18.4091 

(17.795) 

23.3734 

(22.471) 

24.6116 

(23.943) 


45° 

3.1570 

(3.1569) 

8.1435 

(8.0728) 

10.1270 

(10.062) 

18.3843 

(17.79) 

22.9126 

(22.117) 

24.0566 

(23.651) 


60° 

3.0370 

(3.0366) 

8.5855 

(8.4814) 

9.67198 

(8.5911) 

18.3089 

(17.730) 

22.3670 

(21.684) 

23.3533 

(23.160) 


* 2x2, 9-node mesh 

**3x3, 9-node mesh 
t 4x4, 9-node mesh 
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Stricklin et al [9]. The difference between the solutions is mostly in the regions 
of local minimum and maximum. 


6. Two-Layer Cross-Ply Cylindrical Shell Under Uniform Load A cylindrical 
shell with a = b = 5 , R = 10**, h = 0.1" is simply-supported on the four edges. The 
deep shell is laminated by 2 layers (0°/90°) and exerted by a uniform step load 

^ A r 

P = (a P/E^h ) Figure 7 contains a plot of the center deflection versus time for 2- 

D and 3-D elements. The time step used is 6t = 0.1 x 10 - ^ sec. The solutions 
obtained 2-D shell element and 3-D degenerate element are in good agreement. 

7. Four-Layer Angle Ply (45 0 /-45°/45°/-4 5°) Cylindrical Shell Uhder Uniform 
Here we present result for a cylindrical shell which has the same geometry as 

a 

in Problem 3 above. The shell is subjected to a uniform step load P = 50. Figure 
4.8 contains a plot of the center deflection versus time for 2-D and 3-D elements. 
The two elements yield solutions that agree very well in the beginning of the cycle, 
and the 2-D element gives negative values of the deflection at the end of the 
cycle. The discrepency is due to the fact that the 2-D element does not account for 
geometric changes from one time step to next. 

8. Two-Layer Angle-Ply (45°/-45°) Spherical Shell Uhder Uniform Loading 
Consider a spherical shell with a = b = 10", R = 20" and h = 0.1", simply supported 
at four edges and is exerted by a uniform step load. The shell consists of two 
layers, (45°/— 45°). Figure 9 shows the center deflection versus time for P = 50 

A 

and P = 500 with time step 0.2 x lO - ^ sec. For the small load the curve is 
relatively smooth compared to that of the larger load. This is due to the fact that 

A A 

the geometric nonlinearity exhibited at P = 50 is smaller compared to that at P = 
500. 


The results of Problems 3, 4, 6, 7, and 8 should serve as references for future 
investigations. For additional results the reader is referred to Inference 10. 

This completes the discussion of the results. 

CONCLUSIONS 


The present 3-D degenerated element has computational simplicity over a fully 
three-dimensional element and the element accounts for full geometric nonlinearities 
in contrast to 2-D elements based on shell theories. As demonstrated via numerical 
examples, the deflections obtained by the 2-D shell element deviate from those 
obtained by the 3-D element for deep shells. Further, the 3-D element can be used 
to model general shells that are not necessarily doubly-curved. For example, the 
vibration of twisted plates cannot be studied using the 2-D shell element discussed 
in [7]. Of course, the 3-D degenerated element is computationally more demanding 
than the 2-D shell theory element for a given problem. In summary, the present 3-D 
element is an efficient element for the analysis of layered composite plates and 
shells undergoing large displacements and transient motion. 

The 3-D element presented herein can be modified to include thermal stress 
analysis capability and material nonlinearities. While the inclusion of thermal 
stresses is a simple exercise, the inclusion of nonlinear material effects is a 
difficult task. An acceptable material model should be a generalization of Ramberg- 
Osgood relation to a layered anisotropic medium. Another area that requires further 
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Center deflection, wxlO 



Figure 8. Center deflection versus time for four- Figure 9. Center deflection versus time for two- 

layer angle-ply [45 0 /-45°/45 0 /-45°] layer angle-ply [45°/-45°] spherical 

cylindrical shell under uniform load shell under uniform load 


study is the inclusion of damping effects, which are more significant than the shear 
deformation effects. 
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NONLINEAR FINITE ELEMENT ANALYSIS OF SHELLS WITH 

LARGE ASPECT RATIO* 

T. Y. Chang and K. Sawamiphakdi 
The University of Akron 

SUMMARY 


A higher order 'degenerated' shell element with 9-nodes was selected for large 
deformation and post-buckling analysis of thick or thin shells. Elastic-plastic ma- 
terial properties may also be included. A description on the post-buckling analysis 
algorithm is given. Using a square plate, it was demonstrated that the 9-node ele- 
ment does not have shear locking effect even if its aspect ratio was increased to 
the order 10 8 . Two sample problems are given to show the analysis capability of the 
shell element. 


INTRODUCTION 


Research work in finite element analysis applied to plates and shells has en- 
dured for more than twenty years due to the complexity of the problems involved. 
Continuing progress is still being made on topics relating to nonlinear shell analy- 
sis. In particular, recent interest has been focused on large deformation and post- 
buckling behavior of shells. Applications of such analysis problems can be found 
in turbine blades, nuclear vessels and offshore tubular members, etc. 

Development of a finite element procedure for plate or shell analysis can be 
achieved by two distinct approaches: i) using classical shell theories or ii) de- 
riving finite element equations directly from the three dimensional continuum theory. 
Although there are several shell theories of different approximations [1] which are 
useful for linear analysis, they cannot be readily extended to nonlinear cases with 
sufficient generality. Consequently, most of the recent nonlinear shell research 
was concentrated in the latter approach. 

Based on the three-dimensional continuum theory, several different directions 
can be pursued to formulate a shell element. One approach is to deduce a shell ele- 
ment from a 3/D isoparametric solid by imposing necessary kinematic assumptions in 
connection with the small dimension of the shell thickness. Adoption of isoparame- 
tric formulation offers two immediate advantages: i) the requirement of rigid body 

modes is satisfied, and ii) element properties are invariant with reference coordi- 
nates. Several variations of isoparametric-base shell elements have appeared in the 
literature. One class of elements is the so-called 'degenerated' shell family with 
4 to 16 nodes, which was originally proposed by Ahmad, Iron and Zienkiewicz [2] for 
linear shells. Although these elements are quite versatile for extension to nonline- 
ar analysis, several numerical difficulties were experienced. The most notorious 
problem is that the lower order elements exhibit shear-locking phenomenon as the 
thickness of the shell becomes small (or large aspect ratio-element size vs. thick- 
ness). One way to circumvent this problem is to adopt a reduced integration tech- 
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nique for evaluation of element stiffness. An alternate solution is to use higher 
order elements, such as 9- or 16-node Lagrange element. 

In this paper, some of the recent nonlinear analysis results for a 9-node 'de- 
generated' shell element are reported. Nonlinearities considered in our work include 
large deformations, post-buckling behavior and elastic-plastic materials. 


DEGENERATED SHELL ELEMENT 


Detailed description of this element can be found in [3,4] and therefore will 
not be repeated herein. We will only briefly outline this element for the sake of 
completeness. The geometry of the element is circumscribed by its middle surface 
which consists of 9-nodes as shown in Fig. 1. Each node has five degrees of free- 
dom, three translations in the direction of global axes and two rotations about a 
local system. Displacement patterns in the surface of the shell are represented by 
quadratic polynomials. Whereas in the thickness direction, displacements are approx- 
imated by Mind! in's plate assumptions. If the center node of the element is removed, 
it reduces to an 8-node serendipity element. Otherwise, the element is called a 9- 
node Lagrange element. 


SOLUTION METHOD 


The nonlinear shell equations are solved by an incremental tangent stiffness 
approach. For each load increment, either the full or modified Newton-Raphson algo- 
rithm can be optioned in conjunction with secant accelerated iterations. For shells 
exhibiting softening behavior, the modified Newton-Raphson with or without acceler- 
ated iterations was found most effective. On the other hand, for shells exhibiting 
stiffening effect or near instability, the full Newton-Raphson algorithm is necessary 
for obtaining convergent solutions. However, if one is to follow the structural re- 
sponse of a shell beyond its instability point (post-buckling behavior), any of the 
aforementioned algorithms fails to apply due to the singularity of tangent stiffness 
matrix. For this purpose, a different algorithm must be employed. 

There are at least four different methods available for post-buckling analysis: 
i) Artificial spring, ii) specified displacements at nodes, iii) use of current 
stiffness, and iv) constrained arc length. A comprehensive review of these methods 
was given by Ramm [5] and Riks [6]. Of all the methods that have been applied to 
post-buckling analysis of shells, the constrained arc length is most effective due 
to its generality. Actually, in concept this method is equivalent to a displacement 
control analysis, in which numerical instability of a system is circumvented by 
specified boundary displacements. 

There are several ways of defining constrained arc length [7-9,4], but the most 
general definition is as follows. For the i-th iteration of a load increment, we 
calculate an arc length ds by 

ds 2 = a {q i+1 U. +1 )} T {q i+1 (A. +1 )} + 6A? +1 {AR} T {AR} 
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( 1 ) 


= a {q i (X.)} T {q i (X 1 )} + BA* {AR} T {AR} 

and ds must be kept constant, where ds is the arc length at the beginning of the load 
step. In Eq. (1), the following definitions are given: 

q 1 = Incremental displacement vector from time t to t+At after i-th iteration. 

A.. = A load factor after i-th iteration. 

a, 8 = Scaling factors, 0 _< a,8 > 1 . 

For a given problem, the analysis will proceed incrementally with the standard load 
control and the determinant ratio of tangent stiffness is monitored. When the de- 
terminant ratio reaches a small value, i.e. | det KT/det Ko| £ to!., the structure is 
considered to be near unstable. Then the analysis procedure is switched to a con- 
strained arc length method. Thus, the post-buckling behavior of a shell structure 
can be traced without encountering any numerical difficulty. 


ASPECT RATIO 


It is known that the use of 'degenerated' shell elements for thin plate or thin 
shell analysis may give unsatisfactory results due to the so-called shear locking 
phenomenon. This phenomenon was demonstrated for 4- and 8-node elements. One way 
to alleviate this problem is to use a reduced integration scheme [10,11]. However, 
this approach can, at best, postpone the problem and it breaks down when the shell 
thickness is further reduced. An alternate solution to the shear locking problem is 
to adopt higher order Lagrange elements with 9 or 16 nodes. From our study, we found 
that the 9-node element gives very satisfactory results. For discussion purpose, we 
define. 

Aspect Ratio R, = (Largest Element Dimension)/Thickness 

a 

To determine how thin a shell can be modeled by the 9-node elements, a clamped plate 
subjected to uniform load was analyzed by varying the aspect ratio ranging from 10 
to 10 8 . The following cases were considered: 

Case 1. 8-node elements with 2x2 integration order 

Case 2. 8-node elements with 3x3 integration order 

Case 3. 9-node elements with 2x2 integration order 

Case 4. 9-node elements with 3x3 integration order 

Using symmetry condition, one quarter of the plate was sufficiently modeled by a 
4x4 mesh. The plate was loaded well into the large deformation range and the re- 
sults are compared at a load factor (qa 4 /Eh 4 ) = 200. The finite element results for 
all four cases in conjunction with an exact solution are shown in Fig. 2. It is 
clearly seen that the 9-node Lagrange element does not show any shear locking up to 
a ridiculus value of aspect ratio 10®. The problem was also analyzed for a simply- 
supported condition and the same results were obtained. 
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NUMERICAL EXAMPLES 


Two sample problems are presented herein to demonstrate the analysis capability 
of the 9-node Lagrange element together with the post-buckling algorithm described 
in the previous section. 

1. Large Deflection of an Elastic-Plastic Sandwich Cap 

A sandwich spherical cap, shown in Fig. 3, was made of two identical aluminum 
face sheets and a honeycomb core. The face sheets were assumed to be bilinear elas- 
tic-plastic, whereas the core was elastic. The cap was subjected to pressure with 
two variations, i) pressure with constant direction, and ii) pressure always normal 
to the deformed surface (follower pressure). This problem was previously analyzed 
by Sharif i and Popov [12] using two-dimensional axisymmetric elements and the experi- 
mental results were obtained by Lin and Popov [13]. The pressure load was applied 
incrementally in fifteen steps up to 0.8 Pu, where Pu = ultimate pressure of the cap. 
Then the load increment was reduced in half to complete the analysis. As the pres- 
sure was approaching to the ultimate value, the load control analysis was switched 

to constrained displacement method. Throughout the analysis, the modified Newton- 
Raphson iterations with secant acceleration was exercised. It is noted that the se- 
cant acceleration for iteration is activated only if the number of iterations re- 
quired is greater than two (2). On the average, 4-5 iterations per load step were 
used when the cap was becoming structurally unstable. For the case of constant-di- 
rection pressure, our calculated ultimate pressure was found to be 30.15 psi, which 
is fairly close to the buckling load 30.6 psi predicted by Plantema [14] and 29.3 psi 
reported in [12]. For the follower pressure, our calculated ultimate pressure is 
somewhat lower than the constant direction case, i.e. Pu = 27.6 psi. This value is 

compared favorably with the experimental results 27. psi in [13]. The plastic hinge 

was found at a location of 0.8 a from the center of the cap (a = half span of the 
spherical cap). This location is identical to that given in [12]. 

2. Post-Buckling of a Spherical Shell 

A spherical shallow shell, subjected to a concentrated force at the apex and 
supported by fixed hinges, was considered. Two different cases of material proper- 
ties were included: i) linearly elastic material, ii) elastic-perfectly plastic ma- 

terial. This problem was previously analyzed by Argyris et al . [15] using triangular 
shell elements and Parisch [16] using 4-node 'degenerated' shell elements. 

In the case of elastic material, a load increment Ap = 0.1 Pu, Pu = ultimate 
load, was imposed to the shell. Since this structure exhibits prolonged softening 
behavior, the use of modified Newton-Raphson iterations gave considerable difficulty 
in obtaining convergent solution, the full Newton algorithm with 3-5 iterations had 
to be employed. When the load was increased to about 0.9 Pu, the constrained arc 
length method was exercised. From our analysis, the ultimate load was found to be 
Pu = 53.5 lb., slightly higher than the values 51.5 lb. in [16]. 

For elastic-plastic material, our analysis was conducted in four stages in ac- 
cordance with the structural behavior: 

1. Stable region with small deformation (0A in Fig. 4). In this region,, nonline- 
arity of shell is merely caused by elastic-plastic deformation and no numerical dif- 
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ficulty was experienced. With the use of full -Newton-Raphson, only two iterations 
were needed per load step. 

2. Materially unstable, pre-buckling stage (AB). The structure became mildly un- 
stable due to the progression of plastic zone. After the plasticity had spread to 
some extent, internal plastic unloading was taking place as a result of geometric 
change of the shell surface (i.e. large deformation effect). Therefore, there was 

an actual stiffening phenomenon shown in Fig. 4. In this loading range, the determi- 
nant ratio of the shell stiffness was less than 0.1, and consequently convergence 
was difficult to obtain. Analysis was conducted by using the constrained arc length 
method. 

3. Post-buckling stage (BC). The structure is highly unstable and the correspond- 
ing structural stiffness became negative, but still definite. The constrained arc 
length method together with full Newton-Raphson algorithm must be used in order to 
obtain convergent solution. The external load was gradually reduced until a mini- 
mum value was reached (at C). The number of iterations per load step required for 
this region is about 4-5. It is noted that the curvature of the shell surface, as 

a result of large deformation, was changed from convex to concave shape. 

4. Membrane action (CD). After the shell was completely turned upside down, only 
membrane action was present. In this case, the structure has much greater stiffness, 
and hence resumed its stable condition. Correspondingly, the analysis was switched 
back to load control with either modified or full Newton-Raphson iterations. 

The load vs deflection responses for both elastic and elastic-plastic materials 
are compared with those by Argyris [15] and Parisch [16]. Our results correlate 
quite closely with Parisch's solutions, but differ substantially with Argyris 1 2 * * solu- 
tion, especially the post-buckling response of elastic-plastic analysis. 


CONCLUSION 


Applications of a 9-node Lagrange shell element to large deformation and post- 
buckling analysis of elastic-plastic shells have been demonstrated in this paper. 
Based on our study, this element does not exhibit any shear locking effect for a 
square plate problem even though the aspect ratio of the element was increased to 
as much as 10 8 . Obviously, this value is way beyond the range of any real-world 
shell structures. Nevertheless, several additional topics need further attentions. 
Those include the study of its applications to dynamic and creep buckling analysis, 
more efficient ways of evaluating element stiffness to reduce computational effort. 
The latter point may be improved by using a symbolic mathematical manipulation pack- 
age to perform close-form integrations. 
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9- node element with 2x2 integration order 
9-node element with 3 x 3 integration order 


Fig. 2 Max. Deflection of a Square Plate vs. Different Aspect Ratios 







SELF-ADAPTIVE SOLUTION STRATEGIES 
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ABSTRACT 

The paper briefly overviews progress on the development of enhancements to cur- 
rent generation nonlinear finite element algorithms of the incremental Newton-Raphson 
type. The main thrust of the work is to introduce work on new alternative formula- 
tions which lead to improved algorithms which avoid the need for global level up- 
dating and inversion. To quantify the enhanced Newton-Raphson scheme and the new 
alternative algorithm, the results of several benchmarks will be presented. 


INTRODUCTION 

The main thrust of this work is to overview progress on the development and mod- 
ification of algorithms which improve the efficiency and stability of solutions to 
nonlinear finite element (FE) simulations (ref. 1). The emphasis will be threefold 
namely, 1) To review briefly progress on the development of enhancements to current 
generation algorithms of the incremental Newton-Raphson (I NR) type; 2) To introduce 
progress on the development of new alternative algorithmic schemes and; 3) To present 
the results of several benchmark problems. 


OVERVIEW OF PROGRESS ON ENHANCED INR TYPE ALGORITHMS 

While numerous algorithms can be employed to solve the nonlinear algebraic 
equations arising from FE formulations, the most popular scheme is the NR type. 
This follows from the fact that it can handle kinematic and material nonlinearity 
say as in elastic-plastic-creep behavior. Regardless of its adaptability, the 
straight or modified INR suffer from several drawbacks namely: 

1. Cannot handle turning points; 

2. No direct control on successive iterates; 

3. Difficult to ascertain zones of convergence as solution proceeds; 

4. Requires global level updating and assembly; 

5. Global level inversion/pseudo updating required; 

6. Out of core blocking of solution awkward; 

7. Control of individual degree of freedom excursions difficult; 

8. Successive iterations occur on a global level hence control of individual 
degrees of freedom difficult. 

To circumvent the more difficult of the foregoing shortcomings, constraint 
bounds were employed to: 

*Supported by NASA Lewis Under Grant NAG3-54 



1. Handle turning points; 

2. Control successive iterates; 

3. Introduce self-adaptive aspects to algorithm; 

4. Expand regions of stable convergence properties and; 

5. Facilitate evaluation of zones of convergence as solution proceeds. 

In order to facilitate such properties the INR can be reinterpreted as the intersec- 
tion of a tangent extrapolation of the solution curve and a constraint surface. For 
procedures such as the standard NR and those advocated by Wempner (ref. 2) and Riks 
(ref. 3), the constraint surfaces are open hence leading to potential nonintersection. 
In contrast, the use of closed constraint surfaces generally guarantees the intersec- 
tion with the solution curve. While the circular constraint of Crisfield 
(ref. 4) satisfies the closedness criterion, it requires the use of limited step 
size increments. This follows from the manner in which the deflection and force ex- 
cursions are interrelated in normed constraint space. In contrast, the elliptic and 
hyperbolic constraints developed by Padovan and Tovichakchaikul (refs. 5-7) provide 
for more flexible control of successive excursions. 


The formal convergence properties of the constrained methology developed by 
Padovan and Tovichakchaikul has been considered in a very recent series of papers 
(refs. 8-10). This work illustrates various aspects of the formal properties of 
constrained INR schemes namely (refs. 8-10): 

i) The existence of global sized safety zones wherein convergence can be 
guaranteed; 

ii) The use of the various properties of the safety zones to establish self- 
adaptive attributes; 

iii) The existence of quadratic-superl inear convergence rates. 


Based on this work, the elliptically and hyperbolically constrained INR algo- 
rithm have been applied to a wide variety of benchmark problems. Specifically, 

i) Elastic pre-postbuckling behavior of numerous structural configurations 
(ref. 5); 

ii) Elastic-plastic-creep problems including kinematic nonlinearity and poten- 
tial pre-postbuckling (ref. 5-7); 

iii) Nonl inear heat conduction (ref. 11, 12) and; 

iv) Rolling contact problems (ref. 13). 


The results of this work (ref. 5-13) clearly illustrate the enhanced operating 
characteristics of constrained type INR schemes. In particular, the scheme: 


1. Is inherently stable; 

2. Yields significantly improved operating efficiency; 

3. Can be applied to a wide variety of problem types; 

4. Is completely compatible with current architectures of general purpose (GP) 
codes; 

5. Handles turning points; 
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6. Controls successive iterates; 

7. Introduces self-adaptive aspects to INR type algorithms; 

8. Greatly expands regions of stable convergence, and; 

9. Facilitates evaluation of zones of convergence as solution proceeds. 


ALTERNATIVE FORMULATIONS 

The constrained INR approach greatly extends the capabilities to solve nonlinear 
problems. Regardless of this, the NR root of such a formulation still introduces 
difficulties namely 

a) Global updating and assembly; 

b) Global level inversion or pseudo-updating; 

c) Awkward blocking and I/O for out of core problems; 

d) Control of individual degree of freedom excursions 
is difficult; and 

e) Iteration occurs on a global level hence contributing to 
difficulty of controlling individual degrees. 

In the context of the foregoing, an alternative formulation will be sought. The 
main motivations of this work are 

1. Solution scheme should have hierarchial application levels 
(degree of freedom, nodal, elemental, material group, substructural , 
global); 

2. Bypass need of global level inversion; 

3. Bypass need of global level updating; 

4. Develop algorithmic structure enabling simplified 
I/O data flow for out of core problems; 

5. Possess self adaptive attributes enabling hierarchial control 
(locally, globally) and; 

6. Provide continuous hierarchial updating. 

To establish an algorithm which satisfies the preceding requirements, an alter- 
native starting point will be employed. Note, the Galerkin's and virtual work for- 
mulations typically yield algorithms requiring matrix inversion. The use of energy 
type formulations presents difficulties due to it'simplicit form and can fluctuate 
in definiteness during anomalous iterative processes. In contrast, Rayleigh Ritz 
type expressions are positive definite^explicit in form and yield derivatives with a 
simplified structure. 

Based on the standard virtual work principle, the governing FE field equations 
associated with large deformation (small strain large rotation) theory take the 
form (ref. 1). 

R(Y) = / [B*(Y)] T S(Y) dv = F (1) 

V 
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where F is the external load, $ the 2nd Piola-Kirchhoff stress tensor, and Y the 
nodal deflections. Employing (1), the standard Rayleigh-Ritz (least square! ex- 
pression takes the following form 

E <I> - II MY) - F H 2 (2) 

where ||( ) || defines the usual Euclidean norm. To approximate (2), we can employ 
Taylor series. Furthermore, since (2) is positive definite, it has essentially 
quadratic properties about the minimizing solution Y M . In this context we can recast 
(2) as (ref. 14) ~ M 


Taylor approximation; 


eI (Y) = || R(Y m ) - F + [Ky(Y H )](Y - y || 2 

( 3 ) 

Quadratic approximation; 


E "<I) - Y„ + (V - y T [r M ]-’ (Y - Y h ) 

(4) 


such that y M defines the minimum value of E(Y) and [Kj(Y m )] is the tangent stiffness. 

To establish a hierarchial algorithm wherein degree of freedom level updating is 
employed, E* (Y) the quadratic approximation is recast in the form 


E n, <I> ■ 

(R,(]!l> - F, + (v - Y.) T Ky^V,)) 2 + 9, e“ (Y) 

where 


(5) 


E? <P ' Y, + (Y - Y i ) T [r i r 1 (Y - Y,) (6) 

such that is the (i)^ column of [Ky] and 0^ is a scaling factor which controls 
the amount of history of the (i)th degree of freedom admitted per any iteration. 

Employing (5), the localized degree of freedom level algorithm is obtained by 
requiring that the Taylor and quadratic approximation satisfy higher order continuity 
requirements namely 


E H1 <Ii> ■ E H1 <Ii> 


( 7 ) 


d 

W 



<Ii»=57< E ?tl 


(^0) 


(8) 


_df 

dY 2 


(£01 

u i+l 


(Y-)) 




(Y.)) 


(9) 


After extensive manipulations (7-9) yield the following expressions, that is 
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( 10 ) 


Ii+1 ~ li ‘ " F i^ r i^ !Sri 

[r i+i ] = I ([ r i^ - 'i , 1 Cr i ](K T1 ) T (K Ti )[r i 3 T ) (n) 

where is a scaling parameter controlling the modification of [r-j ] . 

The efficient architecture of (10 and 11) follows directly from the fact that 
partitioned matrix multiplication is employed to effect all necessary updating and 
iterating. For example, noting (10) it follows that only a narrow band of [r,-] and 
Kji is employed to update Y^-,. Specifically 

W 5n * 

i- b , 

i 

i+b. 

i-b^ i i+b i 
Similarly noting (11) we see that 
^Ti^ ^Ti ^ = 
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BENCHMARKING 


To benchmark the new algorithm, two highly nonlinear problems have been chosen 
to ascertain the operating characteristics. Figure 1 illustrates the geometry of a 
centrally loaded spherical cap modelled with 8 node quadralateral isoparametric axi- 
symmetric elements. Figure 1 also illustrates the load deflection behavior solved 
via both the standard INR and (10,11). As can be seen, while the standard INR re- 
quired 70 load increments and a total of some 233 iterations, (10,11) achieved a 
converged solution in one load increment involving 7 total cycles through all the 
degrees of freedom. 

As a further demonstration of the use of (10,11), we consider the centrally 
loaded box truss structure given in Figure 2. This structure exhibits pre-post- 
buckling behavior. By imposing linear constraints on the choice of F, the entire 
pre-postbuckling behavior is traced. This example illustrates the ability of (10,11) 
to handle both pre-postbuckling situations wherein transitions in definiteness are 
encountered. 


CONCLUSIONS AND FUTURE OBJECTIVES: ALTERNATIVE FORMULATION 

Based on the foregoing development and benchmarking it follows that the alter- 
native formulation defined by (10 and 11) is: 

1. Stable; 

2. Significantly improves solution efficiency (problem dependent); 

3. Provides for hierarchial application levels; 

4. Eliminates need of inverse; 

5. Architecture of data flow stream-lined; 

6. Can handle problems with turning points (through use of constraints); 
and; 

7. Can handle general material properties. 

The future objectives of work on the alternative algorithm will be several - 
fold namely: 

1. Benchmark algorithm extensively; 

2. Investigate such items as; 

i) Hierarchial application levels; 

ii) Use of localized/globalized constraints to enhance 
solution efficiency; 

iii) Consider use with history dependent media and; 

3. Establish convergence, stability characteristics on formal basis. 
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1 &(W 0 /h) 2 


FIG. 1 Load Deflection Characteristics of Spherical Cap 
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FIG. 2 Load Deflection Characteristics of Box Truss Structure 




ELEMENT-BY-ELEMENT SOLUTION PROCEDURES 


FOR NONLINEAR STRUCTURAL ANALYSIS* 

Thomas J.R. Hughes, James Winget and Itzhak Levit 
Stanford University 


SUMMARY 


Element-by-element approximate factorization procedures are proposed for solv- 
ing the large finite element equation systems which arise in nonlinear structural 
mechanics. Architectural and data base advantages of the present algorithms over 
traditional direct elimination schemes are noted. Results of calculations suggest 
considerable potential for the methods described. 


INTRODUCTION 


Despite the attainment of significant increases in computer storage and speed 
in recent years, many contemporary problems of engineering interest are simply too 
complex to be solved with existing numerical algorithms and presently-available 
hardware. In this paper we address the subject of solving the matrix equations 
arising from finite element spatial discretizations. In Appendix I a brief sketch 
is given of how finite element equation systems emanate from various classes of 
structural mechanics problems. The matrix equations, though sparsely populated, 
still entail enormous storage demands, especially in three-dimensional cases. This 
is the major drawback to matrix-based ("implicit") finite element procedures. The 
types of methods we have developed to deal with this circumvent the need to form 
and factorize large global arrays. These methods have their origins in procedures 
which pervade the numerical analysis literature. Basically, the idea is to replace 
a large, complicated array by a product of simpler arrays. The original concepts 
apparently emanate from the so-called "alternating direction (ADI) methods" (ref’s. 
6, 7, 10, 35). There is a large Russian literature on methods of this type which 
is summarized in the books of Marchuk and Yanenko (ref’s. 31 and 38, resp.) . In 
these works the terminologies used are the "method of fractional steps", the "split- 
ting-up method", and the "method of weak approximation". In the field of computa- 
tional aerodynamics these techniques are often described as "approximate factoriza- 
tion" methods (see e.g. Warming and Beam, ref. 37). The preceding references deal 
primarily with finite difference methods in which the splitting is usually performed 
by decomposing a multi-dimensional partial differential operator into one-dimension- 
al operators. This, of course, places geometrical and topological limitations on 
the discretizations. Generally these methods are used most effectively in the con- 
text of rectangular domains, or domains which are at least topologically equivalent 
to rectangles. When circumstances like this prevail, very large problems can be 
efficiently solved. Progress has been made in developing analogous finite element 


Work performed 


under Grant No. NAG 3-319 


65 



procedures (see ref s. 1, 5, 8, 9, 12—15) . However, these procedures do not retain 
the full geometric and topological versatility of finite element discretizations. 

The methods advocated herein derive many features from the preceding references. 
The approximate factorization aspect of the present approach is facilitated by what 
we feel are the most simple and natural constituents of the finite element process—— 
the individual element arrays . No more than one element array needs to be formed 
and stored at one time and calculations proceed in an element -by-element (EBE) fa- 
shion. There is no geometric or topological restriction imposed by the method, and 
at the same time a remarkably concise computational architecture is achieved. It is 
pointed out herein that the present approach has significant advantages when impli- 
cit-explicit finite element mesh partitions are employed, and, what appears to be 
most significant for the future, the method is amenable to parallel calculations on 
multi-processor computers. 

The idea of element-by-element factorizations was first proposed in Hughes, 

Levit and Winget (ref. 22) in which a transient algorithm for heat conduction was 
developed. Based upon this work, Ortiz, Pinsky and Taylor (ref. 34) constructed a 
novel time-stepping scheme for dynamics. However, our research revealed stringent 
accuracy requirements in certain circumstances, and we were led to reformulate the 
procedure as an iterative linear equation solver (see Hughes, Levit and Winget, ref. 
23) . In this way the usual accuracy and stability properties of standard finite 
element algorithms is attained. The problems that we have applied these procedures 
to are all time dependent and mostly nonlinear. Nour-Omid and Parlett (ref. 33) 
have applied similar procedures to static structures problems and also report en- 
couraging results. 


ITERATIVE ALGORITHMS 


A variety of algorithms may be employed in conjunction with approximately— fac- 
torized arrays. The following has been used in the numerical work presented herein. 

Parabolic Regularization 

The derivation of this algorithm is based upon replacing the algebraic problem 


A x = b 


by a first-order ordinary differential equation whose asymptotic solution is x . 

The terminology "parabolic regularization" is used since the algebraic problem~is 
replaced by what amounts to a spatially-discrete parabolic problem. The ordinary 
differential equation is discretized by backward differences and the implicit opera- 
tor is approximately factorized. Quasi— Newton updates and line searches are employed 
to accelerate convergence. The flowchart summarizes the procedure for symmetric po- 
sitive-def inite arrays. Further details may be found in reference 23. 
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Flowchart of the parabolic regularization (PR) algorithm with line 
search and BFGS updates 


Step 1. Initialization: 

m = 0 , x Q = 0 , r Q = b 

~k = ~k = ~ (loop: k = 1 , 2 



Step 2. 


Step 3 . 


Line search: 


s = Ax t r m /Ax T A Ax 


x . = x + s Ax 

~m+l ~m ~ 


Convergence check: 

lhm+1 11 < 6 


Yes: Return 

No : Continue 


n BFGS^ 


Step 4. Relabel old BFGS vectors: 

~k-l ~ ~k ’ §k-l = Sk ' (loop: k 2 > 3 »•••, n gFGS^ 

Step 5. Calculate new BFGS vectors: 

~ n BFGS = ~m^ ^5 

~ n BFGS = Hm+1 “ ^ " s 2 ^Hm 
Step 6. New search direction: 


5"-5+(fkE> § k <1°°P* k = n BFGS ’ "BFGS " 1 X) 

z B“1 z 

5 <- £ + (g£ ~^~k (loop: k = 1 > 2 »•••» n BFGS ) 

Ax = z 

Step 7. m +■ m + 1 , go to Step 2. 
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The notation in the flowchart is given as follows: m is the iteration counter 

the f^’s and g k ’s are the BFGS vectors; n BFGS is the maximum number of BFGS 
vectors allowed; B is a matrix which approximates A , but is more easily factor- 
ized; s is the search parameter; x m is the mth approximation of x ; r m = b - 
A x m is the corresponding residual; ||r m || is its Euclidean length; and 6 is~a 
preassigned error tolerance. 

Remark 1 . The search parameter in step 2 is determined by minimizing the potential 
energy 


P(s) = - (x + s Ax) T (b - 7 r A(x + s Ax)) (1) 

~m ~ ~ 2. ~ ~m ~ 

Remark 2 . If we ignore the line search and BFGS update ingredients of the PR algo- 
rithm, then the classical Jacobi method is obtained when B is taken to be the 
diagonal of A . 

Remark 3 . Our recent research has indicated that the preconditioned conjugate gra- 
dients method attains faster convergence than the PR algorithm (see ref. 28). In 
addition, conjugate gradients requires only a fixed number of vectors which makes it 
computationally more attractive than the PR algorithm with BFGS updates, because a 
considerable number of BFGS vectors typically need to be stored. For these reasons, 
conjugate gradients is currently our preferred procedure. 


APPROXIMATE FACTORIZATION 


The convergence rate of the algorithm presented in the preceding section 
depends heavily upon the approximating matrix B . It may be noted that if B = A 
then the exact solution x is obtained immediately. Numerous choices for B are 
possible. This subject is explored in reference 28. We limit the following dis- 
cussion to the methods used in computing the results presented herein. To describe 
the procedures employed, we first consider matrices, A , written in the following 
form: 


A = W^I + e A)W^ (2) 

where I is the identity matrix, W is a positive-definite diagonal matrix, e 
is a scalar, and A is a matrix which has the same sparsity pattern as A . A is 
to be thought of as an approximation of A . Specific choices of W , z~ and~ A 
are considered later. The second and final stage of the approximation is to define 

\, 

B = W 2 C W 2 (3) 

where C is an approximation of I + £ A . Various choices are considered below: 
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Two-component Splitting 


Let A be decomposed as follows: 


A ~ A^ + A 2 


(A) 


Then a possible definition of C is 


C = (I + e 


A ) (I + e A 2 ) = 


I + e A + e 2 A A 2 


= I + £ A + 


0(e 2 ) 


(5) 


The last part suggests the nature of the approximation. Computational simplicity is 
gained if and A 2 are very sparse and more easily factorized than A . 


One-pass Multi-component Splitting 


Consider a multi-component sum decomposition of A : 

n 

A = A, (6) 

a ^ 


Let 

n _ 

c = n (I + e A.) 
i=l ~ 

= (I + e A x ) (I + e A 2 ) . . . (I + e A n ) 

= I + £ A + 0(e 2 ) (7) 

Clearly, this is just a straightforward generalization of the two-component splitting. 


Two-pass Multi-component Splitting 


This generalization of the preceding case has qualitative advantages under cer- 
tain circumstances (Marchuk, ref . 31) . Let 


C 


n _ 1 _ 

n (i + f a ) n (i + | a ± ) 

i=l ~ ~ 1 i=n 


(cont'd.) 
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(I +f AjHi +|a 2 ) 


<I + fin> x 


x « + f 6n )( J + f Vl> ••• <1+1*!) 

- Ittlt 0(e 2 ) (8) 

If each Ai is symmetric and positive semi-defiriite, then C is symmetric and 
positive-definite. 


Element-by-element (EBE) Approximate Factorizations 


The EBE approximate factorization is simply a multi-component splitting in 
which the components are the finite element arrays themselves. That is we assume 

n e£ 

A = Ea 6 (9) 

i=l 

where A° is the e th element contribution to A . Then C may be defined by 
either the one-pass or two-pass formulae, viz. 

n e£ 

c = n (I + e A e ) (10) 

e=l 


C = n (I+f?) n (I + I a 6 ) ("Marchuk EBE") 
' 1=1 ' e=n e£ ' 


Remark 1 . We wish to use the term element in the generic sense of a "subdomain 
model", where an element could be an individual finite element or a subassembly of 
elements. Thus we allow limited assembly. Various equivalent terminologies have 
been used to define this concept, such as "substructures" and "superelements" . Sub- 
domain finite element models inherit the symmetry and definiteness properties of the 
global array. Consequently, the remark made after (8) applies. 

Remark 2 . Note that storage demands are vastly less in the EBE case. Only one 
element at a time need be stored and processed. Whether or not it is desirable to 
save factorized element arrays depends upon the availability of high speed RAM, and 
the trade-off between CPU and disk I/O costs. 

Remark 3 . If elements are segregated into non-contiguous subgroups then calculations 
are "parallelizable" . For example, brick-like domains can be decomposed into eight 
non-contiguous element groups (see fig. 1). Because the elements in each subgroup 
have no common degrees-of-freedom, they can be processed in parallel. The eight 
groups, however, need to be processed sequentially. For analogous two-dimensional 
domains, four element groups need to be employed. This feature has significant 
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computational implications for vector and multi-processor machines. 


Remark 4 . It has been our computational experience that if A is symmetric and 
positive-definite* then qualitatively faithful approximate factorizations which 
preserve these properties, such as (11), perform much better than those that do not. 
In the numerical examples presented herein we have employed (11) exclusively. 

Remark 5 . In Reference 28 we have explored the use of reordering the element factors 
and of approximately factorizing the individual element arrays. Preliminary results 
indicate that significant improvements over the basic scheme of (11) can be made. 


SELECTION OF W , e AND A 


In the computations presented herein the following definitions were employed: 

W = diagonal entries of A - (12) 

A = i W - ' 5 A W"* (13) 

Thus 

A = W + A (14) 

This choice is motivated by the derivation of the PR algorithm (see ref. 23). 

Remark 1 . A study performed by Nour-Omid and Parlett (ref. 33) indicates that the 
following procedure, suggested by Winget, is superior: 

A = ^ W* (A - W)W _Js (15) 

where W is given by (12) . This leads to 

A = A (16) 

Remark 2 . The implicit-explicit element concept (see ref f s. 20, 21, 24-27) has a very 
simple and clean implementation within EBE approximate factorizations. Recall that 
an explicit element contributes only its diagonal mass matrix to the coefficient ma- 
trix A . Thus W , according to the preceding definition (i.e., eq. _(12)), totally 
accounts for the explicit element contributions and the corresponding A e? s are 
identically zero. What this means is that explicit elements may be simply omitted 
from the formula for C . In nonlinear problems this opens the way to time-adaptive 
implicit-explicit element partitions. In calculating the element contributions to 
the residual (i.e. "b") a check can be made whether or not the critical time step is 
exceeded for the element. If it is not exceeded, a flag is set to indicate that ele- 
ment contributions to C may be simply ignored. The potential savings in nonlinear 
transient analysis procedures incorporating these ideas is clearly considerable. 
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SAMPLE PROBLEMS 


The computed results were obtained on a VAX computer using single precision 
(32 bits per floating point word) . 


Elastic Cantilever Beam 


The configuration analyzed is shown in figure 2. It represents one-half of a 
plane strain beam modelled with 32 bilinear quadrilateral elements. A lumped mass 
matrix was employed. The loading and boundary conditions are set in accord with an 
exact, static linear elasticity solution (see pp. 35-39, ref. 36). However, here 
the problem is forced dynamically. The beam is assumed initially at rest and all 
loads are applied instantaneously at t « 0 + . In formulating the problem, the 
Newmark algorithm is employed with 3 = % and Y = h (see Appendix I) . With these 
parameters, unconditional stability is attained and no algorithmic damping is intro- 
duced (see ref’s. 11, 17, 21). 

The numerical solution is dominated by response in the fundamental mode. This 
is illustrated in figure 3. At a time step of At = 2.5 x 10~* , an essentially 

exact solution is obtained. At a larger step of At = 2.5 x 10~3 9 a very crude 

approximation of the response is obtained. It is interesting to relate the sizes of 
these s teps to th e critical time step for explicit integration, At cr i t = h m in/ c D = 
hmin/ / (X + 2y) / p = 1.336 x 10~~* 9 and the approximate period of the fundamental 

mode, Tj = .0122 (see table 1). As may be seen, both time steps are far outside the 

range of explicit integration. The larger time step resolves the fundamental mode 
with only 5 steps, and thus is larger than the maximum feasible for this problem. 

In comparing the results of the various methods it is important to keep in mind 
that all methods give identical solutions .* Consequently, the primary basis of com- 
parison is the number of iterations needed to attain the solution. It was found 
that the number of iterations per time step did not vary significantly from one time 
step to another for a given method and specific step size. Results for the first 
time step are presented in table 2. The following observations may be made: In 

general the element-by-element results are superior to Jacobi. Use of line search 
and BFGS updates accelerate convergence. The best results are attained by the ele- 
ment-by-element procedure with line search and BFGS updates. 

It is somewhat surprising that methods (v) and (vi) converge faster at the 
larger time step than at the smaller. At this point we have no explanation for 
this phenomenon. 


Elastic and Elastic-Perf ectly Plastic Cantilever Beam 


The geometrical definition of this problem is identical to the previous one 
except that the entire beam is discretized by a 64 element mesh (the lower part of 
the beam was added to the mesh of figure 2). The boundary conditions were changed 
to the following. 

The convergence tolerance, 6 in step 3 of the flowchart, was taken to be .0l||b||. 
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1 ^( 0 , x 2 , t) = u 2 (0, 0, t) = 0 


,(L, x 2 , t) = Q - (~) ) 


>x 2 € [-c, +c] , t e [0, T] 


Q = 1,000 , T = 0.04 , L = 16 , c = 2 


The boundary tractions are zero on the remaining boundary segments. The tensile 
uniaxial yield stress was taken to be 3,000. Small deformations were assumed and 
the elastic stiffness matrix was used on the left-hand side throughout. The radial- 
return algorithm of reference 30 was employed to integrate the elastic-plastic con- 
stitutive equation. 

Figures 4 and 5 compare the elastic and plastic stress distributions at t = 
.036. A fully developed plastic hinge is present at the root of the beam in the 
plastic case. The elastic critical time step of this problem is At cr it = 1.336 x 
10“5 and the time step used was At = 2.5 x 10“3 = 187.1 • The ^BE method 

with BFGS updates and line searches was employed. The average number of iterations 
for both the elastic and plastic calculations was 4. 


Elastic and Elastic-Perf ectly Plastic Cantilever Beam with a Circular Hole 


The geometric definition of the problem is shown in figure 6. The beam was 
discretized using a 500 element mesh. The following kinematic and stress boundary 
conditions were employed: 


u 1 (0, x 0 , t) = u 9 (0, 0, t) 


,(L, x 2 , t) = q(f )(l - (t) ) 


x 2 6 [-c, +c] , t S [0, T] 


Q = 250 , T = 0.09 , L = 28 , c = 4 


Where not specified to be nonzero, the tractions are zero. The uniaxial yield 
stress was taken to be 1000. A critical time step of At cr £ t = 5.41 x 10~6 was 
calculated on the basis of the smallest element edge length. (The critical time 
step based upon the shortest distance between opposite element edges, which may 
be a more meaningful distance, is less than half this number.) Two time steps 
were employed in the calculations: At = 10 x At cr j_j- and At = 50 x At cr j_t . 

The EBE algorithm with BFGS updates was also employed for this problem. Results 
for the smaller time step converged in 1 iteration, whereas for the latter, 7 iter- 
ations were required on average. 

Figures 7 and 8 show the stress distributions at time t = 0.09 for the elastic 
and plastic solutions. A fully developed plastic hinge is present at the end of the 


73 



beam and a secondary plastic hinge has partially developed in the stress concentra- 
tion zone (plastic regions are shown dashed in fig. 8c). 

We wish to remark upon the contour-line routine used to obtain the results 
presented in this section. The finite element analyzer calculates the stresses at 
the integration points. The data is then extrapolated to the nodal points by means 
of a weighted average of all the integration points in the interior domains of all 
elements connected to the nodal point. The weights are taken as the inverse of the 
distance between the integration point and the nodal point. This type of data 
smoothing ensures that the values obtained at the nodal points will be bounded by 
the data calculated at the integration points which is an essential property when 
plastic stresses are present and ensures continuity of stress contours between ele- 
ment domains. However, this method has some drawbacks. For example, data which is 
anti-symmetric about the neutral axis of the beam, such as <7^ f will result in 
linear distribution of contour lines in all elements having the center line as part 
of their boundary even in the case where these elements have a uniform stress dis- 
tribution (part of a plastic hinge) . Data with symmetry with respect to the neutral 
axis, such as the von Mises stress, does not suffer from this type of smoothing. 


CONCLUDING REMARKS 


In this paper EBE approximate factorization techniques have been proposed and 
compared on test problems. The PR algorithm with BFGS updates performed well, how- 
ever, the need to store BFGS vectors is considered a significant disadvantage, and 
thus the fixed storage requirements and improved convergence characteristics of the 
conjugate gradients algorithm renders it preferable (ref. 28). Improved behavior 
has also been attained in reference 28 by employing reorderings of the element fac- 
tors. 


It should be kept in mind that the EBE concept has been explored herein as a 
finite element, linear equation solving procedure. Although initial attempts at 
directly using EBE ideas to develop time stepping algorithms had some deficiencies 
it may still ultimately prove profitable to couple EBE concepts with the time-step- 
ping loop and even the nonlinear iterative loop. It is interesting to note that the 
inultigrid method found its initial success as a linear equation solver, but in the 
most recent and successful variants the multigrid philosophy permeates all aspects 
of the methodology (see Brandt, ref. 4). 

A step has been taken in the development of EBE solution of finite element equa- 
tion systems. A considerable potential exists for the technique, but much research 
st m remains to be done to bring the methods to fruition. 


APPENDIX I - DERIVATION OF LINEAR ALGEBRAIC SYSTEMS IN THE FINITE ELEMENT 
ANALYSIS OF NONLINEAR MECHANICS PROBLEMS 


Semi-discrete Equations of Nonlinear Mechanics 


Consider the following semi-discrete system 
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M a = F 


(I.D 


where M , a and F represent the (generalized) mass matrix, acceleration vector 
and force vector, respectively. Equation (1.1) may be thought of as arising from a 
finite element discretization of a solid, fluid, structure or combined system. In 
general, M , a and F each depend on time (t) . Explicit characterization of 
M , a and F may be given for particular systems under consideration. 


Nonlinear Structural and Solid Mechanics 


In nonlinear structural and solid mechanics the Lagrangian kinematical descrip- 
tion is frequently adopted. In this case the important kinematical quantities are 
d , the material-particle displacement from a reference configuration; v = d , 
the particle velocity; and a = v - d , the particle acceleration. Dots indicate 
the Lagrangian time-derivative in which the material particle is held fixed. The 
forces are assumed to take the form 


F = F 


ext 


- N 


( 1 . 2 ) 


ext 

where F is the vector of given external forces and N denotes the vector of 

internal forces, which may depend upon d , d and their histories . To make the 
dependence precise, one need introduce equations which define the constitutive (i.e. 
stress-deformation) behavior of the materials in question. These equations vary 
widely in type and complexity. For example, they may be algebraic equations, dif- 
ferential equations or integro-dif f erential equations. In addition, inequality 
constraints may be present, such as in plasticity theory. 


Time Discretization 


To solve the semi-discrete problem, a time-discretization algorithm needs to be 
introduced. For purpose of illustration we shall employ the Newmark family of meth- 
ods (ref. 32). Generalization to other time integrators, such as the Hilber-Hughes- 
Taylor algorithm (ref.’s. 16, 18, 19, 21) which possesses improved properties, may 
be easily facilitated without essential alteration to the following formulation. 

The Newmark "predictors’* are given by 


in + l 

c 

T) ? 
II 

+ At 

V n + 2 (1 - 26)a n 

(1.3) 

~n+l 

= V 
-n 

+ At 

a - Y>5 n 

(1.4) 


where subscripts refer to the step number; At is the time step; d n , y n and a n 
are the approximations to d(t n ) , d(t n ) and d(t n ) , respectively; and 3 and 
Y are parameters which govern the accuracy and stability of the method (ref’s. 11, 
17, 29). 
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Calculations commence with the given initial data (i.e., d Q and v 0 ) and a 0 
which may be calculated from 


”2o 


= F 


ext 


- N, 


( 1 . 5 ) 


if H is diagonal, as is common in structural dynamics, the solution of (1.5) is 
rendered trivial. Otherwise, a factorization, forward reduction and back substitu- 
tion are necessary to obtain aQ . 

In the sequel we are only interested in members of the Newmark family for which 

e > o . 


In each time step a nonlinear algebraic problem arises which may be solved by 
Newton-Raphson and quasi-Newton-type iterative procedures. There are several ways 
of going about this. In the following implementation an algebraic problem is formu- 
lated in which acceleration increments are the unknowns. This form of the implemen- 
tation is useful in that disparate field theories, such as fluid mechanics and heat 
conduction, may be formally considered as special cases. 


i = 0 


Acceleration Formulation 

(i is the iteration counter) 


d (i J 

~n+l 


v (±) 

~n+l 

a (i) 

~n+l 


= d 


n+1 


~n+l 


r (predictor phase) 


= 0 


~ ~ ~n+l ~ ~n+l ~ ~n+l ~n+l (residual, or out-of-balance, force) 

M + yAt + BAt^ (effective mass) 


M Aa = R 


v 


(i+1) 

= a (i) 

•n+1 

~n+l 

.(i+1) 

(i) 

n+1 

~n+l 

(i+1) 

n+1 

= d (i) 
~n+l 


+ Aa 


+ yAt Aa 


+ 3At z Aa 


A 


> (corrector phase) 




( 1 . 6 ) 

( 1 . 7 ) 

( 1 . 8 ) 

( 1 . 9 ) 

( 1 . 10 ) 

(1. 11) 

( 1 . 12 ) 

( 1 . 13 ) 

( 1 . 14 ) 

( 1 . 15 ) 
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If additional iterations are to be performed, i is replaced by i+1 , and 
calculations resume with (1.10). Either a fixed number of iterations may be per- 
formed, or iterating may be terminated when Aa and/or R satisfy preassigned 
convergence conditions. When the iterative phase is completed, the solution at 
step n+1 is defined by the last iterates (viz. d n+ ^ = 4a 11 i Yn+1 = Yn+l 1) > 
and a n +]_ = aO# 1 ) )• At this point, n is replaced by n+1 , and calculations 
for the next time step may begin. 


In practice* d n * v n .and 
along with a^fi^ ; but Vn+l^* 
element level. 


a n are generally saved during the iterative phase 
and d^£j^ may be computed as needed, on the 


The matrices C and K are the tangent damping and tangent stiffness matrices, 
respectively. These are linearized operators associated with N . For example, if 
N is an algebraic function of d and d , then 

K = 9N/9d (1.16) 

and 

C = 9N/9d (1.17) 


Generally in structural and solid mechanics, M , K and C are symmetric, 

M and K are positive-definite, and C is positive semi-definite. 

So-called implicit-explicit mesh partitions (ref's. 2, 3, 20, 21, 24-27) may 
be encompassed by the above formulation simply by excluding explicit element/node 
contributions from the definitions of C and K . A totally explicit formulation 
is attained by ignoring C and K . In these cases it is necessary to employ a 
diagonal mass matrix in explicit regions to attain full computational efficiency . 

It may be observed that the preceding algorithm may be specialized to nonlinear 
statics and linear dynamics and statics: 


Nonlinear Statics 

In this case ignore M and C and set v and a to zero throughout. 

Linear Dynamics 

In this case M , C and K are constant and 

N = C v + K d (1.18) 


Linear Statics 

In this case ignore M and C , set v and a to zero throughout, K is 
constant and 
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N = K d 


(1.19) 


To simplify the writing in the body of this paper we adopt the following nota- 
tions in place of (1.12): 


A x = b (1.20) 

Thus during each step, at each iteration, we wish to solve (1.28) in which A is 
assembled from element arrays, that is 

n e£ 

A = S A e (1.21) 

e=l 
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Table 1 Comparison of time steps used in 
calculations with characteristic 
time scales. 


At 

2.5 x 10 -4 

2.5 x 10 -3 

At 

^ t crit 

18.71 

187.1 

Zl 

At 

48.9 

4.89 


Table 2 Number of iterations required for convergence 
for the problem illustrated in Figure 2. 


At 

Method 

2.5 x 10 -4 

(= 18.71 At 
' cnt y 

2.5 x io ~ 3 
(= 187.1 At crlt ) 

(i) Jacobi 

99 

«■<« 

(ii) Jacobi + LS 

38 

75 

(iii) Jacobi + LS + BFGS 

15 

21 

(iv) EBE 

14 

16 

(v) EBE + LS 

9 

6 

(vi) EBE + LS + BFGS 

5 

4 


Key: LS = line search 

EBE = element-by-element 

Note: (1) No convergence attained after 150 iterations. 
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Figure 1. Decompostion of three-dimensional domain into eight 
groups of brick elements for parallel processing. 
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node A 


Figure 2. Problem definition and finite element mesh 


82 
































At = 2.5X 10"* 



At = 2.5 X 10** 



0,00000 0.00SS3 0.01100 0.01705 0.1 


Figure 3. Vertical displacement of node A 
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b. Plaatic Solution 


b. Plastic Solution 


Figure 4. Stress o x j 


Figure 5. Von Mises Stress 












a. Geometry Definition 
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b. Finite Element Mesh 
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c. Van Mises Stress 
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Figure 7. Stress Distribution 
(elastic solution) 


Figure 8. Stress Distribution 
(plastic solution) 
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ABSTRACT 

The design and implementation of a software system for 
generating finite elements and related computations are 
described. Exact symbolic computational techniques are 
employed to derive strain-displacement matrices and element 
stiffness matrices. Methods for dealing with the exces- 
sive growth of symbolic expressions are discussed. 
Automatic FORTRAN code generation is described with 
emphasis on improving the efficiency of the resultant code. 


1. Introduction 

Recent years have seen increasing interest in using computer-based 
symbolic and algebraic manipulation systems for computations in both 
linear and nonlinear finite element analysis. Application areas 

include the symbolic derivation of stiffness coefficients [ref. 1] , 
[ref. 2], the reduction of tedium in algebraic manipulation, the gen- 
eration of FORTRAN code from symbolic expressions [ref. 3], [ref. 4] 
etc. The benefits and usefulness of such an approach in engine struc- 
ture research are evident. However, several problems need be solved 
before this approach can become widely accepted and practiced: 

(i) the efficiency of the symbolic processor and its ability to handle 
the large expressions associated with practical problems, 

(ii) the interface between a symbolic system and a finite element sys- 
tem on the same computer, and 

(iii) the inefficiencies that are usually associated with automatically 
generated code. 

The MACSYMA system [ref. 5] is a highly sophisticated computer 
system for performing exact symbolic mathematical computations. It 


*Work reported herein has been supported in part by NASA under Grant 
NAG3-298 . 
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provides many high-level commands for symbolic computation such as dif- 
ferentiation, matrix multiplication and matrix inverse. MACSYMA is 
developed at the Massachusetts Institute of Technology and has recently 
been made available on the VAX-11/780 computer which is an affordable 
high-performance minicomputer with a very large address space. The 
version of MACSYMA on the VAX under the UNIX operation system [ref. 6] 
is known as VAXIMA [ref. 71. 


We describe our on-going research on the design and implementation 
of a finite element generator running under the VAXIMA system and our 
experience, so far, in dealing with the problems mentioned earlier. 


2. functional Specifica.ti.Qns and design 


The finite element generator (Generator) as a software system 
should provide a set of predefined functionalities and be constructed 
in such a way as to exhibit prescribed characteristics. These func- 
tional specifications and design goals will then guide the detailed 

program design and implementation of the software system. 

From a functional point of view, the Generator will perform the 
following : 

(1) to assist the user in the symbolic derivation of finite elements; 

(2) to provide routines for a variety of symbolic computations in fin- 
ite element analysis, including linear and non-linear applications 
especially for shells [ref. 8] ; 

(3) to provide easy to use interactive commands for most common opera- 
tions ; 

(4) to allow the mode of operation to range from interactive manual 

control to fully automatic; 

(5) to generate, based on symbolic computations, FORTRAN code in a 
form specified by the user; 

(6) to automatically arrange for generated FORTRAN code to compile, 
link and run with FORTRAN-based finite element analysis packages 
such as the NFAP package [ref. 9] ; 

(7) to provide for easy verification of computational results and 
testing of the code generated. 

In providing the above functions attention must be paid to making 

the system easy to use, modify and extend. Our initial effort is 

focused on the isoparametric element family. Later the system can be 
extended to a wider range of finite elements. 
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2. Generation £f element stiffness matrices 


Symbolic processing can play a very important role in the genera- 
tion of element matrices, especially for higher order elements. As an 
example, we shall describe the automatic processing leading to the 
derivation of the elment stiffness matrix [K] in the isoparametric for- 
mulation. From user supplied input information such as the element 
type, the number of nodes, the nodal degrees of freedom, the displace- 
ment field interpolation polynomial and the material properties matrix 
[Dl , the Generator will derive the shape functions, the strain- 
displacement matrix [B] and the element stiffness matrix [Kl . 

The computation is divided into five logical phases (fig. I) each 
is implemented as a LISP program module running under the VAXIMA sys- 
tem. Aside from certain interface considerations, these modules are 
quite independent and can therefore be implemented and tested 
separately. This allows different people to work on different modules 
at the same time. After the modules are individually tested then they 
can be integrated and verified together. Any problems can be isolated 
to a module and fixed quickly. Detailed description of these phases 
follows. 

2.1. Phase I : define input parameters 

The task of this phase is to interact with the user to define all 
the input names, variables, and values that will be needed later. The 
basic input mode is interactive with the system prompting the user at 
the terminal for needed input information. While the basic input mode 
provides flexibility, the input phase can be tedious. Thus we also 
provide a menu-driven mode where well-known element types together 
with their usual parameter values are pre-defined for user selection. A 
fully user-friendly input phase is a goal of our system. 

The input handling features planned include : 

(1) free format for all input with interactive prompting showing the 
correct input form; 

(2) editing capabilities for correcting typing errors; 

(3) the capability of saving all or part of the input for use later; 

(4) the flexibility of receiving input either interactively or from a 
text file. 

2.2. Phase 11 : Jacobian and (Bl matrix computation 

The strain-displacement matrix [B] is derived from symbolically 
defined shape functions in this phase. Let n be the number of nodes 
then 


® (hj. f h2 , . . . , h n ) 

is the shape function vector whose components are the n shape functions 
h^ through h n . The value for the shape functions will be derived in a 
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later phase. Here we simply compute with the symbolic names. Let r,s 
and t be the natural coordinates in the isoparametric formulation and 
HP be a matrix 


HP 


Hf r 
H, S 
Hf t 


where H,r stands for the partial derivative of H with respect to r. 
The Jacobian J is then 


J = HP . [x,y,z] 

where x stands for the column vector tx,,..., x ] etc. Now the 
inverse, in full symbolic form, of J can be computed as 

1 invj 
det(J) 

By forming the matrix (invj . HP) we can then form the [B] matrix. 


2.2. -Ehasg HI : .s.hape functi o n calculation 

Based on the interpolation polynomials and nodal coordinates the 
shape function vector H is derived and expressed in terms of the 
natural coordinates r.s and t in the isoparametric formulation. Thus 
the explicit values for all h^ and all their partial derivatives with 
respect to r,s and t, needed in HP are computed here. 

l.±. Phase IK: FORTRAN code generation for [B] 

A fortran subroutine for the numerical evaluation of the strain- 
displacement matrix [B] is generated for use with the NFAP package. 
This NFAP package is a large FORTRAN based system for linear and non- 
linear finite element analysis. It is developed and made available to 
us by P. Chang of the University of Akron. It has been modified and 
made to run in FORTRAN 77 under UNIX. 


Assignment statements will be generated to define the components 
of the shape function vector H and the various partial derivatives. 
These variables are then used in assigning values to the FORTRAN array 
corresponding to [B] . FORTRAN code generation is controlled by a "tem- 
plate" file which is a skeleton FORTRAN program with special instruc- 
tions for code generation. 


Details of the FORTRAN code generator will be discussed later. 
2.5,. Pha.ag K : compute and generate FORTRAN code for [K] 


The inverse of the Jacobian, J, appears in [Bl . By keeping the 
inverse of J as INVJ/det(J), the quantity det(J) can be factored from 
[B] and, denoting by [BJ] the matrix [B] thus reduced, we have 


[K] 



(BJ] T . (D) . [BJ] 

det(J) 


dr ds dt 
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The determinant of the Jacobian involves the natural coordinates which 
makes the exact integration in the above formula quite difficult. We 
elect to evaluate det(J) at r=s=t=0 and factor it out of the integral. 
The resulting integrand involves only polynomials in r,s and t which 
are readily integrated. We believe this approximation is reasonable 
and the resulting symbolic expression for [K] will still be more accu- 
rate than results obtainable from numerical quadrature. 

To avoid accumulating large amounts of symbolic expressions, we 
generate the stiffness matrix [K] one entry at a time. Namely, the 
integrand matrix is not formed all at once, instead each entry is com- 
puted and integrated individually. Furthermore, the FORTRAN code of 
each [K] entry is output after it is computed. Usually [K] is sym- 
metric and only the upper triangular part need be computed. Although 
there is a general purpose integration package routine available in 
VAX IMA, we use a specially designed integration program to gain speed 
and efficiency. The integration is organized in such a way as to com- 
bine common subexpressions and produce FORTRAN code which is faster and 
more compact. Details on this later. 


A. The FORTRAN code generator 

A set of LISP programs have been written for generating FORTRAN 
code based on symbolic mathematical expressions derived in VAXIMA. 
This code generator runs under VAXIMA and generates FORTRAN code by 
following a given "template" file (fig. II) . A template file contains 
FORTRAN statements any of which may contain an "active part" which is 
one or several correct VAXIMA statements enclosed in and A 
typical template may look like the following. 


SUBROUTINE STIFF2 ( ST, B, YZ , C , THICK , R, S , ND, KKK ) 
IMPLICIT REAL *8 (A-H,0-Z) 

DIMENSION ST ( 1 ) ,B(4,1) ,YZ(2,1) ,C(4,1) ,SK(18,18) 

C 

C KKK = 1 STIFFNESS CALCULATION 

C KKK = 2 STRAIN-DISPACEMENT MATRIX CALCULATION 

C 

IF ( KKK .EQ. 1 ) GO TO 500 
C Here are some active statements enclosed in U 
{ f77 (det=det j ) , bfortranO } 

RETURN 

500 CONTINUE 

C kfortran generates code for [Kl 
{ kfortran (b,dx, [r,s] ) } 

DO 510 I = 1 , ND 
DO 510 J = 1 , ND 
KL = KL+1 

ST (KL) = SK ( I , J) *THICK 
510 CONTINUE 
RETURN 
END 
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To invoke the FORTRAN generator the following VAXIMA command is used. 
GENFORTRAN (template-file-name, output-file-name) 


GENFORTRAN reads the given template file and passes all characters 
into the output file without modification except for the active parts 
contained in "{" and The VAXIMA commands in the active part will 
be executed by the Generator in the order given and any resultant out- 
put will be directed to the output file as well. Thus, in the above 
example, we can see where the strain-displacement matrix [B] and the 
stiffness matrix EK] are being generated. The output file will be the 
template with all active parts replaced by generated code. This output 
file will therefore be a syntactically correct FORTRAN file ready to 
compile. 

Comment lines are copied into the output file without checking for 
any active parts. This allows the use of { , } inside comment lines. 
Recursive invocation of GENFORTRAN from within a template is allowed. 
This makes it easy to generate a set of related subroutines under the 
control of GENFORTRAN. This arrangement is very flexible and works 
quite well for our purposes. 

When the code generator is fully developed, it will have the abil- 
ity to generate code in separate sessions rather than all at once. 
This means partially generated code can be completed later without re- 
generating what's already been done. When substantial amount of code 
is generated automatically, this feature will become important. To 
help organize the code generated, multiple output files can be automat- 
ically used to contain groupings of subroutines. 

1 . growth and Efficiency of FORTRAN code 

Previous work in using systems such as VAXIMA for finite element 
computation uses user-level programs which does not allow much control 
over the exact manner in which computations are carried out. As a 
result, the ability of handling realistic cases in practice is very 
restricted because of expression growth, a phenomenon in symbolic com- 
putation when intermediate expressions become too large for efficient 
manipulation. 

We use LISP level programs with direct access to internal data 
structures. Thus it is possible to construct programs which will carry 
out computations in such a way as to avoid expression growth as much as 
possible. Therefore, our programs will be able to handle practical 
problems with efficiency. Let us illustrate the control of expression 
growth by the [K] matrix computation. 

First of all, the [B] matrix in the integrand is computed with 
"unevaluated" symbols to keep things small. Thus a typical non-zero 
entry of [B] looks like 


hr 2 hs . y - hr . y hs 2 

where hr 2 = the second component of H, r and y = [y lf . . . ,y n l . But these 
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symbols remain un-evaluated at this point. As stated before, we 
proceed to generate the entries of [K] one at a time to keep expres- 
sions small. In doing this, we apply the formula 

[ki= ££ / ibi i t [D) ij iB) j av 

to collect terms with respect to [D].., Coefficients thus obtained 
are kept in un-expanded form on a lisV which is consulted for dupli- 
cates whenever a new coefficient is generated. This identifies common 
(duplicate) subexpressions in different entries of IK] and keep the 
resulting FORTRAN code compact and more efficient. When a new coeffi- 
cient is formed then it is evaluated and expressed as a polynomial 
involving the natural coordinates r,s and t. Now a special-purpose 
integration routine is used. The integration result for each coeffi- 
cient is converted into FORTRAN code and assigned to a temporary FOR- 
TRAN variable. Thus a typical section of the code for [K] may look 
like the following. 

t30 = -( (16*y3-4*y2-12*yl) *y4+(-12*y2-4*yl) *y3 
1 +4*y2**2+8*yl*y2+4*yl**2)/3.0 

t31 = ( (16*x3-4*x2-12*xl) *y4+(-12*x3+4*x2+8*xl) *y2 
1 +(-4*x3+4*xl) *yl)/3.0 

t32 = ( (16*x4-12*x2-4*xl)*y3+(-4*x4+4*x2) *y2 
1 +(-12*x4+8*x2+4*xl) *yl)/3.0 

t33 = -( (16*x3-4*x2-12*xl) *x4+(-12*x2-4*xl) *x3 
1 +4*x2**2+8*xl*x2+4*xl**2)/3.0 

k (5,7) = 4*(d6*t33+d3*(t32+t31)+dl*t30)/detk 
k (5,8) = 4*(d5*t33+d6*t32+m2*t31+d3*t30)/detk 

In the above, t30 etc. are the temporary variables and dl,d2 etc. are 
entries of the material properties matrix [D] . Without the techniques 
mentioned here, the code for each single IK] entry will require 5 to 8 
continuation cards (for a plane 4-node element) . 

Experiments on the VAX-11/780 with the NFAP package together with 
code for the [ B] and [K] matrices generated show that there is a 10% 
CPU time savings with the above described simplification for the 4-node 
plane element. The savings will be much greater for larger problems. 
Among other things, we are currently studying ways to further simplify 
the expressions for the t's. 
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AUTOMATIC GENERATION OF [B] AND [Kl 
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FORTRAN CODE GENERATION 



Figure II 
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STABILITY AND CONVERGENCE OF UNDERINTEGRATED 
FINITE ELEMENT APPROXIMATIONS* 


J. Tinsley Oden 
The University of Texas 


SUMMARY 

An analysis of the effects of underintegration on the numerical stability and 
convergence characteristics of certain classes of finite element approximations is 
given. Particular attention is given to so-called hourglassing instabilities that 
arise from underintegrating the stiffness matrix entries and checkerboard-instabili- 
ties that arise from underintegrating constraint terms such as those arising from 
incompressibility conditions. A fundamental result reported here is the proof that 
the fully integrated stiffness can be restored in some cases through a post- 
processing operation. 


INTRODUCTION 


In virtually all finite element codes, numerical quadrature is used to evaluate 
integrals defining various matrices in the equations governing the discrete model. 

In recent years, many practitioners underintegrate various entries in the stiffness, 
constraint, or mass matrices so as to either reduce the computational effort required 
to solve a problem or in attempts to improve convergence characteristics of certain 
methods. By "underintegration" we mean the use of quadrature rules of an order lower 
than that required to produce exact integrals of the polynomials appearing in various 
finite element matrices. 

Unfortunately, the use of underintegration frequently leads to serious numerical 
instabilities. Underintegration of the stiffness matrices, for example, may produce 
spurious patterns in the displacement fields which are referred to as "hourglass," 
"chickenwire, " or "keystone" modes, whereas underintegration of constraint terms such 
as those arising from incompressibility conditions may lead to spurious oscillations 
in stress or pressures known as "checkerboard" modes. An example of a checkerboard 
pattern in a stress field in a highly-deformed hyperelastic cylinder under axial com- 
pression, computed using eight-node isoparametric elements and reduced integration of 
an incompressibility constraint, is shown in Fig, 1. The hourglass-instability phe- 
nomena has been discussed by Flannagan and Belytschko (ref. 1) in recent times, where- 
as instabilities due to underintegration of penalized constraints have been studied by 
Oden, Kikuchi, and Song (ref. 2), Oden and Kikuchi (ref. 3), and Oden and Jacquotte 
(refs. 3 and 4) . 

In the present paper, particular attention is given to the issue of hourglass in- 
stabilities in the well-known four-node isoparametric element. 


\his work was supported by NASA Lewis Research Center under Grant NAG3-329. 
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COMPRESSION OF R CYLINDER: Q8N 



CONTOUR SIGMRT 


LORD STEP 1 
LORO FACTOR 0. 2000 


Figure 1 , Computed stress profiles in a compressed incompressible hyperelastic cylinder showing actual 
spurious checkerboard modes. 


STUDY OF A MODEL PROBLEM 


To focus on some key features of hourglass instabilities, it suffices to consider 
a model boundary-value problem defined as follows: 

Find a function u £ H^(^) such that 

a(u,v) = (f,v) for all v in 



where 


a(u, v) 
(f,v) 


| Vu • Vv dxdy 

I fv dxdy 
JQ 


> 


j 


( 2 ) 


Here ft is a smooth bounded domain in the x,y-plane, f is a smooth function de- 
fined on ft, a( • , •) is a continuous bilinear form defined on the usual space Hl(ft) 
of functions with L^-first derivatives, and (•,•) denotes the L -inner product. 
For simplicity, ft will generally be taken as rectangular. Problem (1), of course, 
is a variational statement of the Neumann problem, 

-Au = f in ft . = 0 on 8ft 

; dn 

We next construct a standard, conforming, finite element approximation of (1) on 
a uniform mesh of Q^-elements-(four node elements with bilinear shape functions) , 
the length of the sides of each element being h . The approximate problem is then 

Find u^ 6 H^C H^"(ft) such that 

a ^ U h’ V h^ = ^ f,V h^ f ° r a11 V h in H 

where 

H h = {v h ec°(fi)|v h | fi e Qj/jy} («) 

e 

To evaluate aCu^.v^) in (3), we use a Gaussian quadrature rule of the type. 



a( VV 




( 5 ) 


where 

E = the total number of elements in the mesh 
G = the number of Gauss points (G >1 4) 
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0 

w.. = quadrature weights of element e (1 £ j < 4) 

0 

~j = Gaussian quadrature points for element e (1 £ j £ 4) 

Similarly, (fjV^) may also be computed by Gaussian quadrature — exactly whenever f 
is a polynomial. In this way, we are led to the usual discrete system, 

KU = F ( 6 ) 

where K is the stiffness matrix, U the vector of unknown nodal displacements, and 
F is the load vector. 


Now instead of using (3), we wish to employ an underintegrated bilinear form 


a h ( VV ' £ w e“h ( §e )v h ( !e ) 

e 


(7) 


/N 

with the centroid of element e , corresponding to the use of only one Integra- 

tion point. Then we consider the problem 


Find u h 6 H h 


such that 


VvV 


(f h ,v h > for all v h in 


H 


( 8 ) 


where f h is an approximation of f to be defined later. This leads to the new 
discrete problem 


KU = F (9) 

where F is a modification of F (defined below) designed so as to belong to the 
range of the new underintegrated stiffness matrix K . 

The major questions are: 

1. What is the relationship between U and U (if any)? 

2. Will u^ converge to u as h tends to zero? If not, what can be done 
to enhance the quality of the underintegrated approximation ujj? 

We shall provide an explicit answer to (1) and, while TT^ may not converge, we 
shall show that u^ can be obtained from u^ by a simple post-processing operation! 


HOURGLASS MODES 


We note that (for smooth enough u and v) 


a (u,v) 


Auv dxdy + 

J n 


dn fS vds =< Au > v > 


( 10 ) 
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and ker A = {l } (i.e., the solution to the model problem can be determined only to 
within an arbitrary constant). Likewise, ker K = {(1,1, . . . ,1) }. However, the under- 
integration of a(*,*) enlarges the kernel of A . Indeed, 

ker K = {H,l} (11) 

where 1 = (1,1,..., 1) and 
element^ e , 

K h = 0, 

~e~ 


Thus, the pAeAence o& houJuglaAAtng modeA ahZAeA £ nom the. ^cuZuAe ofi the. mdenJjvtegha- 
ted Att^neAA to adequateZy model the. keAneZ o 6 the opeAatoA defining the gZven boun- 
daAy- value pAoblem. 

Our first result in the direction of studying the hourglass mode is recorded as 
follows : 

Lemma 1 . Let a'(*>*) denote the bilinear f.orm 

a' = a - a h (13) 

where a^,*) is defined in (2) and a^(*,*) is the underintegrated bilinear form 
appearing in (7). Let H £ ker A be a (global) hourglass mode, i.e. 

a^(H,v^) = 0 for all in H* 1 (14) 

Then 

a'(H,w^) = 0 for all w^ in (ker A^)^ - (15) 

where (ker Aj 1 )' L is the orthogonal complement of those functions v, such that 

wv = o for a11 ^ Hh * ^ 

Details of the proof of this and several related results shall be the subject of 
a later paper. It is sufficient to note here that for any v^ we may take the pro- 
jection into ker A^ as 


H is an houAgla&A mode. In particular, for a typical 


K 1 = 0, h 1 = (1, -1, 1, -1), 

^ A/ 


1 = { 1 , 1 , 1 , 1 ) 


( 12 ) 


H 


<H "V OtiT ’ HCkerA h 

Thus, a typical function £ (ker is of the form 

< H ' v h > 


W h ” V h (H,H) 


H 


and one can show from this that (15) holds. 


( 16 ) 
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It is clear from (16) that every v h £ H h is of the form 

+ AH , 6 (ker A^) -1- , A £ JR 

Hence, 

a( vV ■ WV + Xa ' ( v H) 

= *(f,H) + (f,w h ) 

so that 




a ' (Uji’H) 

= (f »H) 

H = hourglass mode 

(17) 



3 ( VV 

= (f,w h ) 

for all £ (ker A^) 

(18) 

Since w, 
h 

can 

be given by (16), 





(f> V ■ 

, f (f *H) 

^ (H,H) 

H, v h ) 


and, therefore, 







WV 

= (f,v h ) - 

(H,H) ^ H,V h^ 






for all v, f H* 1 
h 

(19) 

Hence, in 

(8) we take as f, 
h 

the projection. 




f h - f - ■ 

(f,H) 
(H,H) " 


(20) 

From 

(18), 

we have 






vvv 

- (f,w h ) - 

a ' <u h-V 




VvV 

“ <f h'V 

- a'<vV 





= a h < V T h 

) - a ' (u^,w^) 



Thus, 


100 



V\' i h>V ■ - a ' (u h'V 


(v. ,H) 

= - a ’<vV + TH^r a,(u h’ H) < 21 > 


This result seems to imply that and do not di^QA by an houAgJtcu>&tng 

mode., else the right hand side of (21) would be zero. Thus, the procedures usually 
used to correct underintegrated solutions (by adding multiples of H) may not improve 
solutions. This issue is currently being explored and requires much further study. 


CALCULATION OF STIFFNESS MATRICES 


Confining our attention to a unit element, we introduce the matrices 


T 

s 


{1,-1, -1,1} , s' T 


( 1 , 1 , - 1 ,- 1 } 


h T = {1,-1, 1,-1} , 1 T = {1,1, 1,1} 

T ,123 4, T ,123 4, 

x = {x ,x ,x ,x } , y = {y ,y ,y ,y } 


b[ = f{y 2 -y 4 , y 3 -y r y 4 -y 2 > y r y 3 ) 

,T 1 r \ 

b 2 = 2 {x 4 -x r x 1 -x 3 , x 2 -x 4 , x^x.^ 


\ 


> 


J 


( 22 ) 


where x.,y. = coordinates of node 
i y i 

to a uniform rectangular mesh. If 
element, then the stiffness matrix 
be given by 


i, I s 1,2,3, 4. We no longer confine ourselves 
(£,n) denote coordinates of a unit master Qy- 
K e of a typical finite element can be shown to 


/ Axx T A T Ayy T A T \ 


f / Axx A Ayy A \ 

s.- 8 Nf' + T 1 * 

; \ y Ax y Ax y 

(fi = unit square) 


(23) 


= —(s' 

4 ~ : 


- ss' T ) 


|(sh T -h T s) 


+ ^(hs' T 


s ’h T ) 


(24) 


whereas the underintegrated (one-point) stiffness matrix is 



+ *2# 


(25) 
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and | ft | denotes the area of the element. The rank deficiency of K^ 1 ) follows im- 
mediately from the fact that ~ e 

b^h = 0 , bh = 0 , i = 1,2 (26) 

Oua goal keAe h> to compute exactly K Ioa K ) ^Aom (oa k^) by a 

simple pOAt-pAoceSilng opeAation. That this is poslible follows from thl fact that 


where 


h K h = 16 e 

(1) — T 

K e = Kp ' + e yy 1 

~ C- ~ tl 'V'V 


(27) 

(28) 


and 


Y = h — (h T xb 1 -h T yb_) 

~ ~ Ift I 


(29) 


- 1 
£ 4 


C (£ s T x - ns ,T x) 2 + (£(s T y- n(s' T y)) 2 


ft 


ft + £A + r|B 

o n c 


d£dn 


(30) 


A o = ( >v y 3 } (x r x 2 ) + (y r y 2 )(x 3" X 4 ) 


B o = (y 3 -y 2 ) (x r x 4 } + ty^) (* 2 -x 3 ) 


Note that for parallelograms, z = [(x^xO + (x^-x^) + (yi-y 3 ) + ^ 

whereas for rectangular elements of diameter h e , "e = h^/12|ft^|. “ 4 

The difference matrix, 


K 1 = K -K (1 > 
~e ~e 


(31) 


can be written 


i i5e 

K e= 3—2 


( 32 ) 


Upon assembly, we obtain systems 
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(33) 


KU = F, K^U = F 

r T/N/\ 

F = F - FHH (34) 


in analogy with (3) and (8), respectively. In analogy with (21), 

K (1) (U-U) = | l u e (h^HH - h e /||H||) (35) 

e 

where u^ = h^U/||h e || and H = H/||h||. We can use this result to show that 

||K^ (U-U) || £ ||K f U || (36) 

Thus, the success of the underintegration method will again depend upon how small one 
can make the norm ||k t u|| . This issue is to be the subject of a forthcoming paper. 


CONCLUSIONS 

1) The fully integrated and underintegrated finite element approximations ap- 
parently do not differ by an hourglassing mode. Thus, the introducing of orthogonal 
hourglassing corrections may not be sufficient to produce acceptable solutions from 
underintegrated solutions. This subject requires further study. 

2) On the other hand, it is possible to generate exactly the fully integrated 
stiffness matrix from the underintegrated stiffness matrix using fairly simple post- 
processing operations. An exact formula for such calculations has been presented in 
equation (30). 
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INELASTIC AND DYNAMIC FRACTURE, AND STRESS ANALYSES* 

Satya N. Atluri 

Georgia Institute of Technology 


SUMMARY 


The current work of the author and his students in the areas of (£) large defor- 
mation inelastic stress analysis, (ii) inelastic and dynamic crack propagation, is 
summarized. The salient topics of interest in engine structure analysis that are 
discussed herein include: (a) a new path- independent integral (T)in inelastic frac- 

ture mechanics, (b) analysis of dynamic crack propagation, (c) generalization of 
constitutive relations of inelasticity for finite deformations, (d) complementary 
energy approaches in inelastic analyses, and (e) objectivity of time integration 
schemes in inelastic stress analysis. 


INELASTIC FRACTURE MECHANICS 


The most widely researched topic in elastic-plastic fracture mechanics, in the 
past decade or so, and one that has made certain impressive advances in ductile 
fracture mechanics, has been the now well-known J integral (refs. 1,2). It is also 
well known that J is in fact the component along the crack- line of a vector integral, 
and its significance is in the context of incipient self-similar growth of a crack 
in a (nonlinear) PlcuStic material. In this case, J has the meaning of the rate of 
energy-release per unit of crack extension. In deformation theory of plasticity, 
which precludes unloading and which is mathematically equivalent to a nonlinear 
theory of elasticity, J still characterizes the crack-tip field and is still a path- 
independent integral. However, in this case, J does not have the meaning of energy- 
release rate; it is simply the total potential-energy difference between two identi- 
cal and identically loaded cracked bodies which differ in crack lengths by a differ- 
ential amount. However, in a flow theory of plasticity, even under monotonic load- 
ing, the path- independence of J cannot be theoretically established. Also, under 
aAb'i&UJUiy JLoad hJjstotvLzA , which may include loading and unloading, J is not only not 
path independent, but also it does not have any physical meaning. 

Also, a significant amount of crack growth in a ductile material is necessarily 
accompanied by a significant non-proportional plastic deformation which invalidates 
the deformation- theory of plasticity. Thus, theoretically the validity of J is 
questionable under these circumstances. However, for IJjMSttLd (Xffl OlLYVtk of crack- 
growth, it has been argued (ref. 3) that the strain field undergoes a proportional 
increment due to an increment in applied loading and that J is still a controlling 
parameter. For such situations of J— controlled growth, the concepts of a tearing 
modulus and J-resistance curves have been introduced (ref. 4) to analyze the stabil- 
ity of such growth. Using the aboye concepts, and the related concept of CTOA, en- 
gineering approaches to elastic-plkstic fracture analyses were elaborated upon in 
references 5 and 6. 


Some results of research supported under NASA grant 


NAG3-346. 
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Thus, further research is necessary to understand the mechanics of growth initi- 
ation of cracks in elastic-plastic materials subject to arbitrary loading/unloading 
histories. Also, more-theoretically-valid parameters /criteria are necessary to ana- 
lyze situations wherein there may be substantial amounts of stable crack-growth 
which are no longer within the limits of J-controlled growth. 

Turning to the problem of crack growth in structures operating at elevated tem- 
peratures — the so-called problem of creep crack-growth, numerous experiments have 
recently been undertaken (refs. 7,8 for example) with the purpose of finding a pa- 
rameter which may correlate with the creep crack-propagation rate. Most of these 
investigators considered as candidate parameters, Kj, some form of net section 
stress, or in more recent studies, C*. The parameter C* has been introduced in 
reference 9 and is valid only in pure steady-state creep (e ~ a n ) . This situation 
is characterized by: (i) the body being essentially at steady-state creep condi- 

tions which implies very slow crack propagation and (' ii ) the creep-damage process 
zone being local to and therefore controlled by the crack-tip field. However, sup- 
pose a particular material and geometry result in a crack-propagation rate such that 
elastic strains are not negligible compared to creep rates (ie. , non— steady creep) 
and that at the same time, creep strains are no longer localized to the crack-tip 
region. While C* is not a valid parameter for this case, it appears reasonable to 
expect that crack growth rate is still determined by the local crack-tip field since 
the creep damage process zone is still assumed to be local to the crack-tip. Fur- 
ther, even though C* is a path-independent integral at steady-state creep condi- 
tions, it does not seem to have any physical meaning even under such conditions. 
These and other difficulties with C* have been noted in recent literature (refs. 7, 
8). Thus, there is a need to explore alternatives to C * which remain valid even 
under non-steady creep conditions, which are path- independent integrals and which 
have a well-defined physical meaning so that they may be measured in laboratory 
specimens in an unambiguous manner. 

With the above objectives, the author has studied (ref. 10) certain path-inde- 
pendent integrals, of relevance in the presence of cracks, in nonlinear elastic and 
inelastic solids. The hypothesized material properties included: (£) finite and 

infinitesimal elasticity, (ii) rate- independent incremental flow theory of plastic- 
ity, and (Hi) rate-sensitive behavior including viscoplasticity and creep. In each 
case, finite deformations were considered, along with the effects of body forces, 
material acceleration, and arbitrary traction/displacement conditions on the crack 
front. Also, the physical interpretations of each of the integrals either in terms 
of crack-tip energy release rates or simply the energy-rate differences in two com- 
parison cracked bodies were explored. Several differences between the presently 
obtained results and those considered well established in literature were pointed 
out and discussed. 

We omit the mathematical details of the development but present here the end re- 
sults for incremental, path- independent, vector integrals which are of relevance in 
ductile fracture mechanics under arbitrary load histories. Further, for simplicity, 
we consider: (£) quasi-static crack growth initiation and stable crack growth and 

(i£) the deformations to be small, so that the distinctions between various stress 
measures disappear. Consider a generic increment of external load (load-control) 
from an initial state with stresses t and initial displacements _u. Let the incre- 
mental stresses be At and the incremental displacements be Au.. We consider two 
paths T 1 and ^ surrounding the crack- tip. We allow the incremental potential AU 
for At to be an explicit function of the location of the material particle through 
parameters a (which can be 1 for plastic loading and 0 for elastic unloading) and 
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g (the strain-hardening parameter). Thus, there can be plastic loading/or elastic 
unloading from a plastic state/or pure elastic-loading taking place arbitrarily near 
the crack-tip as well as elsewhere in the body. Then it has been shown (refs. 10,11) 
that the following path- independent integrals prevail: 


(AT) p 


lr> j.r [x:Ae + AU)N - N*x*Ae - N*Ax*e]ds 
Jr i +r c ~ ~ ------- 

+ /„ „ [Vt + iVAr):Ae - (Ve + iVAe):Ax]dv 

t 1 


Jr^+r 


Ae + AU)N - N*x*Ae - N*Ax*e]ds 


2 c 


+ 


and (AT*) 

“ P 


I v [(Vx + ^VAx):Ae - (Ve + zVAe):Ax]dv 
= [(x:Ae + AU)N - N*x*Ae - N*Ax*e]ds 


= J^[(x: Ae 

+ jv-v [A ~ T:( ? 


+ AU)N - N*x*Ae - N*Ax*e]ds 
(Ve + £VAe) - Ae:(Vx + iVAx)]dv 


( 1 ) 


( 2 ) 


In the above V t is the total volume of the cracked body, and V 2 are the volumes 
enclosed by and T e is a path arbitrarily close to the crack- tip, and V £ is 
the volume enclosed by r e . It can be shown (ref, 11) that path- independence in the 
sense of equations (1,2) prevails even when there is arbitrary loading/unloading 
occuring in and V 2 and thus between and V^. In the above, ( )•( ) implies a 
tensor inner product and ( ) : ( ) implies a trace (see ref. 10 for details). It has 
been shown (ref. 11) that in Mode I, (AT^)p has the meaning: (i) it is the incre- 

mental area corresponding to an incremental load between the load-deflection curves 
of two identical and identically loaded cracked bodies with crack lengths differing 
by a differential amount, (ii) this physical meaning remains valid in loading as 
well as unloading, and (Hi) under monotonic loading when near-proportionality pre- 
vails near the crack-tip, SAT = J (the summation is for all increments of loading, 
using rate theory of plasticity). In reference 11, a numerical experiment was also 
discussed wherein a compact tension specimen was loaded to a certain load level, 
unloaded to a compressive load, and reloaded again. The observations are: (i) dur- 

ing loading E(AT^)p = J and (ii) (AT]_)p remains strictly path- independent even for 
subsequent unloading and reloading, whereas, J not only has no meaning during un- 
loading, but its numerical value varies widely (±50%) over different paths. 

Thus, while much further work remains to be done, it is evident that the inte- 
grals in equations (1,2) provide a theoretical basis for elastic-plastic fracture 
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mechanics under arbitrary loading and when the material behavior is characterized by 
realistic rate theories as opposed to a simple deformation theory. Analogous deri- 
vations of path- independent incremental vector integral (AT) C valid in QZneJvxt un- 
steady C/iee.p, and a discussion of its merits as compared to the C* integral, have 
been presented elaborately by Stonesifer and Atluri (refs. 12-14). 


DYNAMIC FRACTURE MECHANICS 


We first consider linear elastodynamic crack propagation in a self-similar fash- 
ion. Consider, for instance, a two-dimensional problem. Let the velocity of crack 
propagation be denoted by the vector C, with modulus C. Consider a small loop r e , 
of radius e, centered around the crack- tip at time t, and assume that the crack-tip 
moves by an amount (dtC) into T e at time t+dt, ie. , e > Cdt without loss of general- 
ity. Let the unit outward normal to r e be N, and let the kinetic energy density be 
T(= ipyv, where v is the absolute velocity of a material particle and p the mass 
density). Note that v is singular near the propagating crack-tip. From studies in 
reference 10, it is seen that the rate of energy release rate G per unit of elasto- 
dynamic crack growth is given by: 

G = (1/C) Lt C'f [ (W+T)N - t-e]dT = (1/C)C*G (3) 

e->o e 


Now we consider the problem of analyzing crack-propagation in an arbitrary body, 
the shape of which and the loading on which, we suppose, preclude any possibility of 
an analytical solution. Suppose that we have to use a numerical solution. Such a 
numerical solution may be based on a "propagating singular element" within which the 
asymptotic mixed mode solution is embedded; and hence the K-factors can be evaluated 
directly, as demonstrated by the author’s group (refs. 15-18). However, in order to 
use a simple numerical procedure, say using distorted (singular) isoparametric fi- 
nite elements or non-singular isoparametric elements, it is convenient to have 
available path- independent integrals which have the same meaning as the energy re- 
lease rate G of equation (3). If so, the integral can be evaluated on a path that 
is far removed from the crack-tip and hence is insensitive to the details of model- 
ing of crack-tip stress-strain fields. 

Such a path-independent integral has been derived based on general conservation 
laws, as shown in references 10 and 19, and is given below: 

•e]dr + f v _ v [p(a-F)*e - VT]dv 

£ 


r-fr 


[ (w+t)n - t 


L 


[ (w + -w' 


■)N + + 


(t + -t' 


)N + - 


t*e]dT = G 


(4) 


where (+) and (-) refer, arbitrarily, to the "upper" and "lower" crack faces; N + = 
-N - is the unit normal to r c + ; and t are prescribed tractions on the crack face. 
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The practical application of the above new path- independent integral has been 
demonstrated in references 20 and 21. 


GENERALIZATION OF CONSTITUTIVE RELATIONS 
TO FINITE DEFORMATIONS 


At a symposium on finite strain plasticity held at Stanford University in 1981, 
Nagtegaal and de Jong (ref. 22) presented some interesting results for stresses 
generated by simple shear of elastic-plastic and rigid-plastic materials which ex- 
hibit anisotropic hardening. In the equation for the rate of change of the shift 
tensor a, they used the Zaremba-Jaumann-Noll rigid body rate of a. They found the 
rather spurious result that the shear stress is oscillatory in time. 

The above 'anamoly' has prompted a series of investigations by Lee and his 
associates (refs. 23,24). As a 'remedy 1 , Lee et al. (refs. 23,24) suggest the use 
of a 'modified' Jaumann derivative of a in the evolution equation for a. While the 
use of such a modified derivative in the specific problem of simple shear has been 
illustrated in references 23 and 24, its generalization to three-dimensional, non- 
homogeneous cases is not yet fully developed. It has been suggested in reference 
23 that M a complete investigation of the micro-mechanics and the structures of poss- 
ible macroscopic constitutive relations will no doubt be needed to fully understand 
this phenomenon and to generate a fully-tested theory". 

Recently the author has shown (ref. 25) that the anamolies as described above 
are not peculiar to the anisotropic plasticity alone; similar behaviour in finite 
shear may result even in the case of hypo-elasticity and classical isotropic-harden- 
ing plasticity theory. Thus, in reference 25, the central problem of 'generalizing' 
to the finite deformation case, of the constitutive relations of infinitesimal 
strain theories of elasticity (hypoelasticity) , and of classical plasticity with 
isotropic or kinematic hardening was discussed in detail. It has been shown in ref- 
erence 25 that the current controversies surrounding the choice of stress rate in 
the finite-strain generalizations of the constitutive relations and the anamolies 
surrounding the kinematic hardening plasticity theory are easily resolvable. 

Here we briefly illustrate the case of kinematic hardening plasticity. In this 
case, considering a simple J 2 theory, the current yield surface can be represented 
by: 


f= (o' - a'): (o' - a') - — - = 0 


(5) 


where O ' denotes the deviator of a tensor, g is the Kirchhoff stress tensor, and 
a is the back stress. The normality condition is: 


e P = X(3f/3o) 


( 6 ) 


We define, as in reference 25, generalized 'objective* rates, o* and a* of q and a 
respectively, as follows: 
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(7) 


2 * = §° ~ Y y(?*e) “ YgCg'S) 

and a* = a 0 - Yy(a*e) - Yg(e*a) (8) 

where a 0 can be any one of the currently-used stress-rates. For instance. 


n m T 

= ? - e*a - a« e (Oldroyd) (9) 

a 0 = a m + + a-e (Cotter-Rivlin) (10) 

a° = a m _ e *g 4- a *e (mixed-1) (11) 

T T 

gO = a m + e *o - o • e (mixed-2) (12) 

= a m + a*u - w*a (Zaremba-Jaumann-Noll) (13) 


In equations (7-13), e is the velocity gradient, e the strain-rate, oj is the spin, 
and d m is the non-objective material rate. Since, for a classical continuum without 
body couples, the material rate <5r m should be symmetric, the generalized stress rates 
a* and a* are defined in reference 25, such that: 

Py = Yy = Yg if is that of Zaremba-Jaumann-Noll (14) 

Py = (Yy + 1) = (Yg+1) if or 0 is that of Oldroyd (15) 

Py = (Yy~l) = (Yg“l) if is that of Cotter-Rivlin (16) 

Py = (Yy“l) = (Yg+1) if that of mixed-1 (17) 

Py = (Yg“l) = (Yy+1) if is that of mixed-2 (18) 

The linear kinematic hardening theory of Prager, as generalized to finite deforma- 
tions, states that: 

.* P 

ot = eg ( 19 ) 
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where c is a constant of proportionality. It can be shown (ref. 25) that the plas- 
tic strain-rate can be written, for the yield function given in equation (5), as: 


e P = —ry [(a'-a'):a*](a'-«') 
2ca z ----- 


( 20 ) 


Using equations (7,8), (14-18), and (20), the material rates d and a can be writ- 
ten as: 


d m = ce 4- wot - a*a) + y^(a*e + e»a) 


( 21 ) 


a m = 2ye + X(e:I)I - 2y _ 9 [e: (o'-a ') ] (a'-a ') 


i - 
l 


3y 


(c+2y)c^ 


+ to *a - a*to + y^(a*e + e *a) 


( 22 ) 


In a finite shear problem, we may prescribe velocities: 


V 1 = 2 “ X 2 V 2 = ° 


(23) 


such that the strain- rate e and the spin to are, respectively: 



"0 to 0" 


0 to 0 

£ = 

— 1 

O O 

O O 

l°i 

and oj = 

-to 0 0 
. 0 0 0. 


Lee et al. (refs. 23,24) consider the case when: (i) Y y and Yg are zero in equa- 

tions (7) and (8); (ii) yy = 0 in equations (14,21,22) ie., they consider the case 
when cr° is the Zaremba-Jaumann-Noll rate. The obtained solutions are oscillatory in 
references (23,24). Thus, the generalization of Lee et al. (ref. 23,24) -Li -impnopzn. 
On the other hand, it has been shown by the author (ref. 25) that if yy is ±1 in equa- 
tions (21,22), the obtalnzd 6ohutiorU> ion. 6tn.Z66QA oJiz non-oAcittatony. Furthermore, 
it is seen that the condition yy = ±1 can be satisfied even when any one of the 
stress-rates as in equations (9-13) are used. The way to do this is to UA Z thz 
gznznatizzd &tA.< 2 A&-natQA a* and a* aA In zquatiom { 7, S ) and to zhooiz y 7 and Yg a 6 
tn zquattoni [14- IS). 


COMPLEMENTARY ENERGY APPROACHES IN INELASTIC ANALYSIS 


Recently Reed and Atluri (refs. 26-31) presented comprehensive studies of a new 
hybrid-stress finite element algorithm, suitable for analyses of large, quasistatic. 
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deformations. The algorithm is based on a new complementary energy prin- 
ciple for rate type materials given earlier by Atluri (ref. 32). The principal var- 
iables in this formulation are the nominal stress-rate and spin, and the resulting 
finite element equations are discrete versions of the equations of compatibility and 
angular momentum balance. 

The algorithm produces true rates, time derivatives, as opposed to ’increments'. 
There results a complete separation of the boundary value. problem (for stress-rate 
and velocity) and the initial value problem (for total stress and deformation); 
hence, their numerical treatments are essentially independent. In references 30-31, 
a comprehensive discussion of the numerical treatment of the boundary value problem 
and detailed examination of the initial value problem covering the topics of effi- 
ciency, stability, and objectivity were presented. Several problems dealing with 
homogeneous as well as inhomogeneous large deformations of inelastic materials were 
solved, and the solutions discussed in detail. In the following, we present a brief 
outline of the methodology. 

The boundary- value problem for the rates of stress and spin is governed by a 
variational principle corresponding to the stationarity of the functional: 

- ix:(o)*m) + £:w}dV + n*t*v dS (25) 

v 

where w is the spin, and t is the rate of the first-Piola stress as referred to the 
current configuration. In the above, R is the rate of complementary energy density, 
t is the current Cauchy stress, v are prescribed velocities at the boundary— segment 
S v of the solid, n is a unit outward normal to S v , and V is the volume of the solid. 
In equation (25), the following definitions and a priori constraints apply: 


V*t + pb = 0 in V 


ort = Vx<j) + t 


(26) 


n* t 



at S 

a 


(27) 


9R 

3r 


e 


(28) 


• • T • T 

r = i(t + T*U! + 0J • T + t) 


(29) 


In the above, V is the gradient operator in the current configuration of the body, 

and t are prescribed body forces; $ is a first-order (once-dif ferentiable) stress- 

function, t b is a particular solution for t; and e is the strain rate. The Euler- 

Lagrange equations that follow from 6n =6 are: 

c 

(e+u)*t + t = t + T*(g-ho ) ( 30 ) 
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( 31 ) 


v = v at S 
- - v 


Vv = 


(e+u)^ 


(32) 


Equation (30) is the angular momentum balance condition, equation (31) is the veloc- 
ity boundary condition, and (32) is the kinematic compatibility condition for the 
rate problem. 

In the finite element formulation, we assume for each element: 


NQ 

v - £ 5 ± q* ; N 

i=l 


= isoparametric shape functions 


(33) 


co = 


NW 

i=l 


where QW. + QW. =0 
~ ~ 1 


(34) 


NT 

t = X/ QT.3^ + t^ where QT. = V x <J> and V*t^ = -pb 
- i-i - 1 * - - 1 - ~ s - 


(35) 


where NQ is the number of velocity parameters, NW is the number of spin parameters, 
and NT is the number of stress parameters. 

Use of equations (33-35) in (25) can be shown to lead (ref. 30) to the algebraic 
system: 


Kv* = f 


(36) 


where v* is the vector of global nodal velocities. Solution of (36) leads to global 
velocities, from which the rate t and the spin co in each element can be computed. 

We are thus led to initial value problem to integrate v, t, and y in time to find 
the solution for displacement, total stress, and rotation at each material point. 

In references (30-31) an elaborate study of various time-stepping schemes for 
integration of the above initial value problem as to their efficiency, stability, 
and objectivity was presented. In the following, we limit ourselves to a brief 
description of an important and interesting aspect of this time integration — which 
pertains to the question of maintaining "objectivity” of numerical integration. 
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OBJECTIVITY OF NUMERICAL INTEGRATION 


In order to be called objective., a numerical approximation for a physical entity 
must transform between two observer frames according to the same rule as the entity 
itself. An algorithm which produces an objective approximation will itself be call- 
ed objective. 

Suppose that as an objective rate one were to consider the Zaremba-Jaumann-Noll 
rate, ie. , 


Do 

~~ + O' *03 — 0) •O' 
~ ~ ~ - 


(37) 


Suppose also that the solution at time t = T is known, ie. , cr(x) = a at t = T. If 
the usual Euler integration approximation were used, ie. , 


a(x+At) = a(x) + (d° - a* a) +aj*a)At 


(38) 


the result will be in error, since (I + a>At) is not an orthogonal tensor. 

Recently two algorithms were presented (refs. 33-34) which preserve the stress- 
integration and give a proper expression that is different from equation (38). In 
reference 30, Reed and Atluri presented an approach that is slightly different from 
that of references 33-34 in that numerical objectivity is treated as a general ques- 
tion without specific reference to stress— integration. This approach involves the 
concept of the so-called Jaumann integral (ref. 35) using which objectivity can be 
preserved in any numerical integration. Here, we present the use of the Jaumann in- 
tegral in stress-integration. 

The Jaumann stress— integral, which can be viewed as operation inverse to the 
Jaumann-differential of equation (37), can be written (ref. 30) as: 

T(t) = J~ 1 (t)9(t)-x(x)-Q T (t) +J [J t (4)Q(t)-9 T (?)*a°(4)-Q(c)*Q T (t)]dc (39) 


where x is the true stress, and Q(t) is the solution of the differential equation: 


9 ( t) =m(t)*g(t) ; Q(x) = I (40) 

Now let tg be the time such that 

t 0 = x + 0At 0 < 6 < 1 ( 41 ) 

we replace equation (40) by an appfwXAjncutZ differential equation: 
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Q(t) - u 0 *g(t) ; w 0 =w(t + 0At) 


(42) 


and treat y 0 as a constant in equation (42). The 2.xa,ct solution of equation (42) 
would be: 


9 0 = 9d 0 ) = ( ex P* ‘i)e t 0^*J 


where 


Now let 


(! " c °SC0 0 t 0 ) 2 

I H oj . + (o. 

~e a). ~0 


0 


0 


w e = 


Q 1 = Q(t 1 ) = Q(x + At) 


(43) 

(44) 

(45) 


Then, an appWXAJncution to the Jaumann integral, equation (39), can be written as: 


t(t + At) = J T 1 (t 1 ) •Q 1 *t*Q^ + [(At)J 1 (t 0 )Q 1 *Q 0 *a°*g 0 *Q^] ; 0 < 0 < 1 


T vi 


(46) 


where t^ = x + At. For 0 = % we obtain the so-called "generalized midpoint rule", 
which agrees with that presented earlier in reference 34. The formula of reference 
33 cannot be recovered since they use a midpoint constitutive evaluation but fail 
to properly rotate the resulting stress increment. The present concept of "objec- 
tive integration" has the advantage that it provides a general framework for any 
time- integration. 
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SUMMARY 


The major objective of this work has been to adapt general nonlinear finite- 
element codes for the purpose of analyzing the dynamics of gas turbine engines. In 
particular, this adaptation required the development of a squeeze-film damper ele- 
ment software package and its implantation into a representative current generation 
code. The ADINA code was selected because of our prior use of it and familiarity 
with its internal structure and logic. This objective has been met and our results 
indicate that such use of general purpose codes is a viable alternative to special- 
ized codes for general dynamics analysis of engines. 


INTRODUCTION 


There is currently a considerable interest and level of activity in developing 
computational schemes to predict general engine dynamic behavior. The general feel- 
ing among researchers working on engine vibration problems is that various modes of 
operation, such as blade-loss events, require a high level of analysis sophistication 
to realistically model the engine. Proper account of system nonlinearities (parti- 
cularly at the bearings, dampers and rubs) appears to be necessary if analytical 
predictions are to be realistic. The approach described in this paper seeks to make 
use of already proven general finite- element nonlinear time-transient computer codes 
which are available on the open market. 

Present-day jet engine configurations have evolved to a substantial degree 
through a trial-and-error process involving extensive testing. There are many fun- 
damental dynamic phenomena which take place within these engines for which basic 
description and understanding have yet to be generated. Nonetheless, they work well. 
Modern aircraft engines are typical of current high-technology products in which the 
recently acquired computing capabilities of today are being used to better under- 
stand and improve what is already designed, built and operating. 

A better understanding of the basic dynamic characteristics of existing and new 
engine configurations is a prerequisite for producing acceptable engine efficiencies 
on advanced configurations (i.e., smaller rotor/stator running clearances). Also, 
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a better definition of engine dynamic response would more than likely provide valu- 
able information and insights leading to reduced maintenance and overhaul costs on 
existing configurations. Furthermore, application of advanced engine dynamic simu- 
lation methods could potentially provide a considerable cost reduction in the deve- 
lopment of new engine configurations by eliminating some of the trial— and— error 
process done with engine hardware development. 

The emergence of advanced finite element codes, such as NASTRAN, NONSAP, MARC, 
ADINA, ANSYS and ABAQUS and related algorithmic advances, have placed comprehensive 
engine system dynamic analyses within reasonable reach. What remains to be done is 
to develop new component element software to properly model engine rotor/stator 
interactive components, such as the squeeze-film damper, within the algorithmic logic 
of already proven finite-element codes. This is the major mission of this work. 

For good reasons, aircraft gas turbine engines use rolling element bearings 
exclusively. This design philosophy has, until recent years, deprived engines of 
beneficial damping inherent in many other types of rotating machinery where fluid- 
film journal bearings are used. The implementation of squeeze-film dampers in 
recent engine designs has now provided engine designers with an effective means of 
vibration energy dissipation. The net result is that engines with squeeze-film 
dampers are less sensitive to residual rotor imbalance and better able to control 
vibration and transmitted force levels resulting from various excitation sources 
within the engine. 

The field of rotor dynamics has evolved to its present state primarily through 
the solution to problems in classes of machinery older than aircraft engines. In 
most other types of rotating machinery (e.g., steam turbines, centrifugal pumps and 
compressors, fans, generators, motors, etc.) the rotor can be adequately modeled as 
an Euler or Timoshenko beam [1]. In addition, the support structure holding each 
bearing can often be adequately modeled as a separate mass-damping-stiffness path to 
ground (i.e., to the inertial frame). Also, for most purposes, bearing dynamic pro- 
perties are characterized as stiffness and damping elements, linearized for small 
vibration amplitudes about some static equilibrium state. With few exceptions (e.g., 
Hibner [2]), it is this level of sophistication that has been utilized for the most 
part in rotor-dynamics analyses of aircraft engines. 

Present day aircraft engines are structurally far more complex than most other 
types of rotating machinery. The multi-shaft configuration, plus the fact that the 
shafts are thin rotating shells, creates unique but significant complicating diffe- 
rences between aircraft engines and other machinery. Also, the stator structural 
support at each rotor bearing represents anything but a separate mass-damper-stif f- 
ness path to an inertial frame. In fact, setting the inertial frame for the engine 
is not a simple matter when the full range of in-service maneuvers is realized. Dy- 
namic paths between different bearings exist not only through the rotor but through 
several other paths within the nonrotating engine structure, i.e., a "multi-level," 
"multi-branch" system. As many as eight significant "levels” have been identified. 

Adams [3] demonstrated the use of a free-free modal vector space for the rotor 
in the numerical integration of the motion equation with nonlinear journal bearings 
and base-motion inputs. The resulting computer code (ROTNL) has been applied to 
large steam turbine blade-out unbalance, instability limit cycles and base-motion 
excitation of nuclear machinery (earthquake and underwater explosive shock). This 
same approach was later applied by Gallardo et al . , [4] in developing the TETRA code 
for double-spool shaft gas turbine engines. The approach used in references [3] and 
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[4] uses a standard code (i.e., NASTRAN) to generate the free-free undamped normal 
modes and natural frequencies for rotor and stator configurations. Nonlinear forces 
and connections between substructures are carried on the right hand side of the 
equations of motion, and updated at each time step of the particular numerical inte- 
gration scheme. In this approach, a special purpose code is used to perform the 
integration of the motion equations, and handle input and output processing. Consi- 
derable flexibility and continual program advancement is therefore convenient. Also, 
this provides easier control of the basic mathematical model. 

On the other hand, the major advantage of the approach summarized in our work is 
that full advantage can be taken of the multiplicity of computational and modeling 
advancements typical of the capability in current generation nonlinear general pur- 
pose finite- element codes. The disadvantage is that the general purpose code is more 
of a ’’black box” than one’s own special purpose code. 


TIME-TRANSIENT NONLINEAR DYNAMIC ANALYSES 


In recent years it has become evident that an important class of engine dynamic 
phenomena can not be studied without accounting for the highly nonlinear forces pro- 
duced at bearings /dampers, labyrinths and other close-running rotor/stator clearances 
under large amplitude vibrations. In such cases, linear theory typically predicts 
vibration amplitudes larger than the actual running clearances. Furthermore, impor- 
tant vibratory phenomena, such as subharmonic resonance and motion limit cycles, are 
"filtered" out of the problem with a linear model, giving grossly erroneous predic- 
tions, qualitatively as well as quantitatively [3]. 

With few exceptions, nonlinear dynamics problems must be solved numerically as 
time-transient responses, whether the sought answer is a steady-state periodic 
motion or is strictly a transient phenomenon. The problem is mathematically catego- 
rized as an initial value problem in which the displacements and velocities of the 
complete system must all be specified at the beginning of the transient. From that 
point forward in time, the equations of motion are numerically integrated (known as 
"marching) as far in time as one wishes to study the system motions and forces. If 
the system is dynamically stable, the transient motion dies oilt yielding the steady 
state response which in a system with a periodic force excitation will be periodic 
motion. In a stable system with no time-varying force excitation, the transient will 
die out as the system comes to rest at one of its stable static equilibrium positions. 
If the system is unstable, the transient does not die out but continues to grow in 
time unless or until some nonlinear mechanism in the system limits the motion to 
what is frequently called a "limit cycle" [3]. 

In order to study the general dynamical characteristics of aircraft engines, 
nonlinear dynamics computational schemes are required. The approach taken is to de- 
velop software packages to model engine components which are not typically found on 
dynamical structures and therefore are not already built into existing nonlinear 
finite element structural dynamics computer codes. The initial effort has concen- 
trated on developing such a software package for squeeze-film bearing dampers. 
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OVERALL APPROACH-INTERACTIVE ELEMENTS 


Considering the typical engine structural complexities, an improved computa- 
tional approach is necessary if a proper transient/ steady-state model is to be deve- 
loped for gas turbine engines. In this approach, it appears that the finite-element 
method is one of the attractive modeling techniques for such problems. Its inherent 
capabilities include features essential to modern engines: 1) automatically handles 
multi-branch, multi-level structures in a more direct and efficient manner than fle- 
xibilities approaches, 2) well-suited to handle nonlinearities associated with struc- 
tural kinematic and kinetic effects [7], 3) easily accommodates various types of 
boundary and constraint conditions, and 4) easily accommodates material nonisotropy 
and nonlinearity [7,8]. A body of established and proven algorithms are available 
which can handle these various important effects [7,9] as well as geometric comple- 
xities (e.g., beam, plate, 2-D and 3-D elements [10]. 

The required features which are presently not available with general purpose 
finite— element codes are provisions to handle rotor/stator interactive forces origi- 
nating from squeeze-film dampers, seals and rub/impact events. Presented herein are 
the results of an effort to develop a squeeze-film damper computer software package 
which can be "plugged" into existing finite-element codes. This work is detailed in 
references [5] and [6]. 


DEVELOPMENT AND IMPLANTATION OF SQUEEZE-FILM DAMPER ELEMENT 


The basic implant approach is quite simple as shown in Figure 1. In addition 
to the standard structural element input provided for the particular general purpose 
code (i.e., host program) being used, an additional small block of fixed inputs is 
provided for each damper in the model. In the computational stream for numerical 
integration of the motion equations, the damper software routine is called at each 
time step of the marching scheme, and provides the total instantaneous interactive 
force generated in the fluid film of the squeeze-film damper. This instantaneous 
force is then applied to the structures as an "external" force. The damper also 
provides an option to compute the instantaneous tangent stiffness and damping matri- 
ces for the damper if such is needed for improvement to the numerical integration 
scheme of the host program. Details of the formulation and solution technique used 
in the squeeze film software package are detailed in reference [5] • 

To "experiment" with the damper element independent of a general purpose code, 
a small four degree-of-freedom host program was written which treated the rotor and 
stator each as point masses with two degrees of freedom apiece, "joined" to each 
other only by the squeeze-film damper (see Figure 2) . Sample runs are shown in 
Figures 3, 4, 5 and 6. 

As detailed in reference [6], the damper element was successfully "implanted" 
into the ADINA code and several case studies were made. These are summarized here 
in Figures 7, 8 and 9 for a multi-bearing simple rotor and different types of 
transient loading. 
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EXPERIMENTAL CONFIRMATION 


A two-bearing two-damper test rig was designed and fabricated. The rotor con- 
sists of a slender constant diameter shaft with a concentrated disk. mast at the 
center of the bearing span (see Figure 10). Extensive testing was performed with 
this set up and rotor vibrational orbits were measured using Bently proximity probes 
and signal conditioning instrumentation. We are currently completing this work by 
making runs for this configuration using the ADINA program with squeeze-film soft- 
ware package implanted. We anticipate completing our correlations of experimental 
and computational results in the immediate future. 
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Figure 1 Input/Output of Damper Code, 
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Figure 2 Simple 2-Mass 4-Degree of Freedom Test Model. 
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(b) Rotor Test Bed With Squeeze- Film Dampers 

Figure 10 Existing Experimental Test Facility Rotor 
Dynamical Studies With Squeeze-Film Dampers. 
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SUMMARY 


An integrated approach to nonlinear analysis of high-temperature composites 
in turbine blade/vane applications is presented. The overall strategy of this 
approach and the key elements comprising this approach are summarized. Preliminary 
results for a tungsten- fiber-reinforced superalloy (TFRS) composite are discussed. 


INTRODUCTION 


A major thrust of current aeropropulsion research is directed at improving 
the performance/efficiency of aircraft gas turbine engines. One area of research 
where efforts have been focused concerns the potential use of advanced and high- 
temperature composite materials for critical engine hot-section applications 
such as high-pressure turbine (HPT) blade/vane components. The temperature- 
limited capability of the conventional materials (mostly nickel- and cobalt- 
based monolithic superalloys) used in HPT blades/vanes is a key factor affecting 
the ability to augment engine performance/efficiency. 

The ability to increase the maximum engine operating cycle temperature 
and, thus, improve engine performance/efficiency would be greatly enhanced 
by improving the elevated temperature capability of turbine blade/vane 
materials. For this reason a particular family of the refractory-fiber 
metal-matrix category of high-temperature composites, namely the tungsten- 
fiber-reinforced superalloys (TFRS) , are being investigated (refs. 1-2) as 
a potential alternative to conventional monolithic superalloys for use in 
HPT blade/vane components.- Moreover, a specific TFRS composite has been 
identified (ref. 3) as having an excellent combination of complementary 
properties to make it a viable candidate as a first-generation composite 
turbine blade/vane material. 

From a qualitative viewpoint, the potential benefit of using high- 
temperature composite materials in engine hot-section applications, such 
as TFRS for HPT blades/vanes, is evident. Despite the promising outlook, 
a dilemma exists regarding the lack of adequate methods available to 
quantitatively assess a high- temperature composite in an engine hot-section 
application. For example, any present attempt to analyze the behavior of 
TFRS in a HPT blade/vane application would be forced to rely on experimental 
evaluations which are economically impractical. A definite need exists 
for an analytical capability with which to make a quantitative assessment 
of the mechanical performance and structural integrity of high -temperature 
composites, such as TFRS, in engine applications, such as HPT blades/vanes. 
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In view of the above discussion, a formal program was initiated at the 
NASA Lewis Research Center to develop a structural/stress analysis capability 
specifically tailored for the nonlinear analysis of TFRS composite HPT blades/ 
vanes. Although the capability has been tailored for a specific application, 
the overall approach of this analytical capability is general and as such is 
applicable to the analysis of practically any composite structure where 
material nonlinearity is a concern. 

The emphasis throughout the development of this structural/stress analysis 
capability for TFRS composite HPT blade/vane components has been to incorporate 
the physics governing this application into the analytical formulation of the 
problem. The blade/vane components of the HPT section operate in severe 
environments which subjects them to complex and cyclic thermomechanical loading 
conditions. These loads give rise to highly nonlinear material behavior which, 
for the case of TFRS composites, is assumed to be attributable to a stress- 
temperature-time dependency of the constituent (fiber/matrix/interphase) 
material properties. Furthermore, in the extended service elevated temperature 
environment of a HPT blade/vane, TFRS composites are subject to a metallurgical 
phenomenon known as fiber degradation which involves a chemical interaction 
at the fiber-matrix interface and results in the generation of an interphase 
of material. 

To develop a nonlinear structural/stress analysis capability for TFRS 
composite HPT blade/vane components, the COBSTRAN computer code was used as 
a foundation. COBSTRAN is a software package developed as an in-house research 
tool at the NASA Lewis Research Center to perform linear-elastic structural 
analyses of composite blade structures. The resulting nonlinear version of 
COBSTRAN incorporates the appropriate nonlinear thermoviscoplastic (TVP) 
material relationships, fiber degradation expression, composite micromechanics 
equations, laminate theory, and global structural analysis to relate the global 
nonlinear structural response of a HPT blade/vane component to the behavior 
of the TFRS composite constituent materials. 

Development of a nonlinear version of COBSTRAN was undertaken as an 
evolutionary task. The intent is to "fine-tune" nonlinear COBSTRAN through 
an iterative process of empirical investigation and analytical methodolgy 
refinement. During this iterative process, nonlinear COBSTRAN will be used 
as a tool to conduct parametric/sensitivity studies to isolate the critical 
factors affecting the behavior of TFRS composites in HPT blade/vane applications. 
Information from the parametric studies will be used to identify appropriate 
experimental investigations to further examine the critical factors that 
warrant accountability in a structural analysis. 


ANALYTICAL APPROACH 


The overall approach taken here for the nonlinear structural analysis of 
a TFRS composite HPT blade/vane component is illustrated in the schematic in 
Figure 1 . This schematic is representative of the general solution strategy 
of nonlinear COBSTRAN. This approach is referred to as an upward- integrated 
top-down-structured analysis. Figure 1 illustrates the levels of heirarchy 
involved in this approach and summarizes the key elements comprising this 
approach . 
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The circular arrangement of Figure 1 is intended to illustrate the 
incremental nature of a nonlinear structural analysis and the iterative 
process involved in the solution of such a problem. The cycle corresponds 
to the analysis of each distinct load increment. The solution process for 
a nonlinear analysis by this approach involves iteration at two levels. 

A primary iteration occurs for the load increments of the total solution 
range. For each increment, a secondary iteration occurs to establish 
equilibrium, in the integrated sense, of the global structural response. 


Nonlinear Thermoviscoplastic Relationships 


As pointed out earlier, material nonlinearity for a TFRS composite 
HPT blade/vane is attributed to a stress-temperature-time dependency of the 
constituent material properties. To account for or model the material non- 
linearity in a structural analysis, the nonlinear TVP relationships shown 
in Figure 2 have been implemented into nonlinear COBSTRAN for the purpose 
of updating the constituent material properties at each iteration during 
the analysis solution process. 

As shown in Figure 2 the nonlinear TVP expressions provide a relationship 
between the current value (P) and reference value (P Q ) of a material property 
according to the current and reference state of the constituent material as 
described by the indicated field variables. The properties being modified 
include; normal moduli (E-q, E 22 ) > average heat capacity (C) , thermal 
conductivities (K^, K 22 ) > thermal expansion coefficients C a n» a 22^ » an ^ 
remaining strengths (S;q, S 22 > S 22 ) • 

These particular nonlinear TVP relationships are presented with the 
realization that they do not represent an exact or complete mathematical 
model for describing the nonlinear material behavior for a composite structure. 
The expressions in Figure 2 represent an attempt to account for the material 
nonlinearity involved in the application of TFRS composites in a HPT blade/ 
vane compenents by treating the problem at the fundamental level. 


Composite Micromechanics 


In short, composite micromechanics is the branch of composite mechanics 
which provides the formal structure to relate the behavior of a composite 
lamina to the behavior of the constituents. Composite micromechanics equations 
have been derived to predict; the material properties (mechanical, thermal, and 
uniaxial strength) of a unidirectional fiber-reinforced lamina based on the 
corresponding properties of the constituent materials, and, the distribution 
of microstresses in the constituent materials resulting from the stress state 
occurring in a lamina. 

In the lamina transverse directions (2- and 3-directions, see Figure 3), 
the distribution of constituent microstresses 0^22 > a 12» a 23^ non-un if orm 
through the thickness of the ply. By virtue of the nonlinear TVP relationships 
discussed earlier, this intralaminar non-uniformity of transverse microstresses 
results in a non-uniformity of transverse material properties for the constituents. 
This intrlaminar non-uniformity is illustrated in Figure 3 in terms of the 
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regions or zones used to characterize the non-uniformity. Figure 3 also 
lists the equations used to calculate average constituent transverse properties. 

The composite micromechanics equations for predicting lamina mechanical 
properties are shown in Figure 4 together with a schematic defining the material 
property coordinate axis system. Similar equations exist for predicting 
lamina thermal properties and uniaxial strengths as well as for predicting 
the distribution of constituent microstresses from the corresponding lamina 
stresses. 

In a sense, the composite micromechanics equations presented here can 
be thought of as a first level of approximation for idealizing an inhomo- 
geneous laminated TFRS composite HPT blade/vane component as a pseudo- 
homogeneous structure. 


Laminate Theory 


In short, laminate theory (laminated plate theory) provides the formal 
structure to describe the behavior of a laminated composite structure in terms 
of the individual laminae comprising the laminate. A cursory summary of the 
governing equations of classical laminated plate theory is given in Figure 5. 
Just as the composite micromechanics equations were considered as a first level 
of approximation, laminate theory can be thought of as a second level of 
approximation in idealizing an inhomogeneous TFRS composite HPT blade/vane 
component as a pseudo-homogeneous structure. 


Global Structural Analysis 


The global finite element formulation of a general linear structural 
analysis problem is summarized by the equations shown in Figure 6. Although 
this collection of equations presents a superficial view of the finite element 
method of structural analysis, these equations are representative of the 
governing equations of structural analysis as implemented in a general purpose 
finite element code such as NASTRAN. 

The basis of the integrated approach presented here for the nonlinear 
analysis of a TFRS composite HPT blade/vane component involves the application 
of the linear finite element formulation summarized in Figure 6 in an incre- 
mental manner. 


PRELIMINARY RESULTS AND DISCUSSION 


Results were obtained from the analysis of a single unidirectional 
TFRS composite lamina using the integrated approach just presented as 
implemented into nonlinear COBSTRAN. Details of the analysis are left 
out here as the purpose of presenting these results is merely to: illus- 
trate the capabilities of this integrated analytical approach by showing 
the types of information obtainable; identify trends of material behavior 
being accounted for and traced in a nonlinear analysis by this integrated 
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approach, and; point out the possible significance of a fabrication 
process in conducting a nonlinear analysis of a composite engine 
structure such as a turbine blade/vane. 

The results in Figure 7 were obtained from the analysis of a 
single unidirectional TFRS composite lamina subjected to a temperature 
load representative of the cooling portion of a typical fabrication 
process. The first plot illustrates the large longitudinal normal 
residual stresses (a-^) existing in the constituents at the end of 
this fabrication process. These stresses are the result of the 
difference in thermal expansion coefficients between the constituents. 
The second plot shows the variation of longitudinal normal modulus 
(E^j) for the constituents and ply during the fabrication process. 

Note that the modulus for the matrix initially increases but then 
decreases continuously for the duration of cooling. This reflects 
the relative influence or effect of the individual dependency terms 
in the nonlinear TVP relationships presented earlier. At first, 
the temperature-dependency dominates and, as might be expected, the 
modulus increases to a point where the residual stress becomes 
large enough such that the stress -dependency becomes dominant. The 
results for the fiber and interphase (degraded zone) can be explained 
in a similar manner. The result for the ply, as determined from the 
corresponding micromechanics equation, reflects the behavior of the 
fiber and matrix as the volume fraction of interphase was nominal in 
this case. 

The third plot in Figure 7 illustrates the variation of remaining 
strength available in the ply and constituents during the fabrication 
process. The results are consistent with the corresponding nonlinear 
TVP relationship for strength. Finally, the last plot shows the 
build-up of transverse normal residual microstress (0 22 ) the 

constituents similar to that shown for a^. The point to be noted 
here is the difference in microstress from region to region for the 
matrix and interphase (degraded zone) . This reflects the intra- 
laminar non -uniformity discussed earlier. 

The conclusion to be drawn from these results is that perhaps 
the fabrication process for individual laminae, laminate panels, and 
a composite component itself should not be neglected. 
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SYMBOLS 


volume fraction 

exponential constants (equations. Figure 2) 
total time 

distance to lamina (equations. Figure 5) 

laminate membrane stiffness matrix 

heat capacity, laminate coupled membrane -bendind 
stiffness matrix, or global structural damping 
matrix 

constants (equations, Figure 2) 

initial fiber diameter and current fiber diameter 
or laminate bending stiffness matrix 

normal modulus or material property matrix 

shear modulus 

thermal conductivity or global structural 
stiffness matrix 

laminate moment resultants or global structural 
mass matrix 

laminate stress resultants 

thermal and mechanical load cycles (equations. Figure 2) 

initial and current value for a property (equations. 
Figure 2) 

lamina transformation matrix 

initial, current remaining, cumulative, and fracture 
strengths (equations. Figure 2) 

initial, current, and melting temperatures (equations. 
Figure 2) 

thermal expansion coefficients 
strain, and laminate midplane strains 
laminate curvatures 


Poisson's ratio 



<v a 


initial, and current stress 


° 0 > o , o. 


H 


0) 


subscripts : 
c 

F, M, D, L 
1, 2, 3 
x, y, z 


initial, current, and reference maximum stress 
rates (equations. Figure 2) 

vibration frequency 


composite or laminate 
fiber, matrix, interphase, and lamina 
lamina material property coordinate axes 
global structural axes 
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COMPONENT 



Composite 

Micromechanics 



Figure 1 - Upward-integrated top-down-structured approach for nonlinear 
analysis of a composite engine structure 
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Figure 2 - Nonlinear thermoviscoplastic (TVP) relationships to characterize 

composite material nonlinearity in terms of stress-temperature-time 
dependent constituent properties 
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AVERAGE TRANSVERSE PROPERTIES (E 22 , G iz , & n , K 22 , ot 22 ) 


matrix: p„= [p M w + (l-|) P M <» + yr F (|) P„ ( ‘>] 


INTERPHASE : P p = 


7J [f- (IWl 


Figure 3 - Various regions used to characterize intralaminar non-uniformity 
of transverse constituent properties and microstresses 
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(F) FIBER 
(M) MATRIX 
(D) INTERPHASE 
(L) PLY 
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Figure 4 - Composite micromechanics equations for predicting average 

mechanical properties of a unidirectional fiber-reinforced ply 




Y 


/ 


PLY CONSTITUTIVE RELATIONSHIP : = [E]^(T} l + aT {a} L 

PLY strain: {e\= [P] L <Je°} c ^ \{^} c )> 

'g PLY stress: {<t} l = [E] L <({e} L +- ATjot}^ 




LAMINATE CONSTITUTIVE RELATIONSHIP-* 



{N t ,M T } c 





UMINATE STIFFNESS: [A, C, 0] c = I [R][ [EjjR^ <( (2 L - i L „), J *’,))> 


Figure 5 - Summary of the governing equations of laminate theory 
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Figure 6 - Summary of the global finite element formulation of a 
general linear structural analysis problem 
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Figure 7 - Preliminary results of an analysis of a single unidirectional lamina 
subjected to the cooling half of an hypothetical fabrication process; 
variation of longitudinal normal stress (ctu) 
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Figure 7 - continued; variation of longitudinal normal modulus (Ejj) 
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Figure 7 - continued; variation of longitudinal tensile strength (S 11T ) 
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THREE-DIMENSIONAL STRESS ANALYSIS USING THE BOUNDARY ELEMENT METHOD* 


Raymond B, Wilson 
Pratt & Whitney Aircraft 

Prasanta K. Banerjee 
State Universuty of New York - Buffalo 


SUMMARY 


The boundary element method will be extended (as part of the NASA Inelastic 
Analysis Methods program) to the three-dimensional stress analysis of gas turbine 
engine hot section components. This paper outlines the analytical basis of the 
method (as developed in elasticity), summarizes its numerical implementation and 
indicates the approaches to be followed in extending the method to include 
inelastic material response. 


NOTATION 


V 

S 



5 ij 

e . 

J 

E 

v 

M 

a 

e 


interior of three-dimensional body 

surface of three-dimensional body 

displacement vector 

strain tensor 

stress tensor 

unit normal vector 

traction vector 

body force distribution 

Kronecker delta 

arbitrary unit vector 

Young's modulus 

Poisson’s ratio 

shear modulus 

coefficient of thermal expansion 
temperature 


*Work done in support of NASA contract NAS3-23697. 
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INTRODUCTION 


It has been established that the durability of hot section components is a 
major driver of gas turbine engine costs, in design, manufacture and maintenance. 
Major savings will be realized if the lives of these components can be extended and 
their cost reduced. The NASA HOST (Hot Section Structural Technology) program is 
supporting a variety of approaches to these problems. As part of the HOST program a 
need has been identified for an improved inelastic three-dimensional stress 
analysis capability which is efficient enough for use in the design cycle. NASA 
Contract NAS3-23697 (3-D Inelastic Analysis Methods for Hot Section Components), 
recently awarded to the Commercial Engineering Division of Pratt & Whitney 
Aircraft, is directed at the development of analysis tools to meet this need. 

The analysis requirements are severe. Hot section components, and the related 
stress analyses, are characterized by: 

1. Extreme geometrical complexity. The components typically include surfaces 
of very general shape. Small structural details must be accurately 
represented, since these are frequently life limiting locations. The 
components usually contain holes of various sizes to allow for air and 
fuel flow and these can both weaken the overall structure and act as local 
stress risers. It may also be necessary to account for the presence of 
single or multiple cracks within the structure. 

2. Loading complexity. Loads typically include gas path and cooling air 
pressures, mechanical loads due to other components, thermal loads 
(derived from unsteady temperature distributions) and, in some cases, 
extremely significant centrifugal loads. Durability analysis for hot 
section components requires consideration of the stress-strain history 
throughout a flight cycle, and thus all of the boundary conditions and 
other loads are potentially time dependent. 

3. Complex material response. In all cases the elastic material response in 
hot section components is inhomogeneous due to the high temperatures 
involved and the temperature dependence of elastic modulus. In some cases 
the materials are anisotropic. Inelastic response usually occurs, at least 
in local regions, over some part of the flight cycle. In some cases global 
inelastic response (plasticity, creep or both) of the component is 
encountered. 

The requirements of the design process cannot readily be met using a single 
analysis tool. The early phases of the design process require rapid and inexpensive 
analyses of a variety of designs. In the final phase of design analysis a very 
accurate determination of the stress-strain state is required for final life 
prediction, but greater analysis time and cost are acceptable. 

Much of the work in the present program will be devoted to the development of 
simplified (beam, shell) analysis methods for use in early design and to the 
development of special purpose finite element modules for the efficient analysis of 
hot section components. Significant effort will also be devoted to the extension of 
the boundary element method to the three-dimensional stress analysis of hot section 
components. This will provide a complementary method which can be used as an 
alternative to finite element analysis. It is expected that, for some problem 
classes, accuracy and efficiency will be improved relative to the use of finite 
element analysis. 
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The remainder of this paper reviews the analytical development of the elastic 
boundary element method, briefly discusses its numerical implementation, gives an 
example of its use for a rather complex structure and indicates the approach to be 
adopted in extending the method to inelastic analysis. 


ELASTIC ANALYTICAL FORMULATION 


The boundary element method is a numerical method applicable to the solution 
of problems in homogeneous elasticity. It is similar to the finite element method 
in utilizing piecewise approximations to geometry and field quantities. These 
approximations allow the effective treatment of complex structures such as those 
encountered in gas turbine engines. The boundary element method is based on an 
entirely different mathematical development than is the finite element method, and 
this leads to profound differences in the numerical implementation. The key feature 
of this development is the use of classical results in mechanics and potential 
theory to reduce the dimensionality of the problem. Thus a three-dimensional 
elasticity problem is solved solely in terms of the surface of the body. The 
mathematical development is discussed below. Further details may be found in 
reference 1. 

The starting point for the development is the Reciprocal Theorem (ref. 2) 



a & 

where ( a ij> e ij) anc * (<Jij> € ij) are t ^ ie stresses and strains derived 

from any two displacement fields satisfying the equilibrium and compatibility 

equations. Use of the strain displacement equations and the divergence theorem 
allows the recasting of the Reciprocal Theorem in an alternate form involving 
surface integrals of the displacement and traction fields, and volume integrals of 
the displacement fields and body force distributions. 


Jt± uj_ dS + ^fjL u-f_ dV = Jt i Uj_ dS + Jf 


av 


( 2 ) 


The displacement states chosen are completely arbitrary, except for the 
requirement that equilibrium be satisfied. For isotropic, homogeneous materials one 
can identify the (*) solution with the displacement field due to a unit force 
applied in the direction ej at the point y^ (the Kelvin solution, ref. 1). 

Formally this is the displacement field satisfying 



axj 


8 ij 


8(x,y) ej = fi 


(3) 
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and it can be represented in the form 


* 


* 


U£ (x) = G.JJ (x,y) ej (y) 


( 4 ) 


where 


G i i = f(3 - 4 v ) 

3 16t r/i (1 - v)r l * 

(r = |x - y| ) 



+ 


(Xi - y^Xj -y.,) 
r2 


(5) 


The strains, stresses and (on the surface, S) the tractions can be represented 
in similar forms. 


6 ij ®ijk e k 

* _ T * 

°i j - T ijk e k 

t i = a ij n j = ^ijk n j e k 88 ^ik e k 


where the tensors ^ijk an ^ ^ik are derived from the displacement 

tensor, G^j, by appropriate differentiation and use of Hooke f s Law. A similar 
approach is available for anisotropic materials, although closed form 
representations are not generally available. 

By substituting u*, t* and f* in the Reciprocal Work Theorem 
and utilizing the formal definition of fj, an identity (Somigliana 1 s identity) 
can be derived which expresses displacements at an arbitrary interior point in 
terms of the displacements and tractions on the surface. Appropriate 
differentiation and use of Hooke's Law allow derivation of similar results for 
stresses and strains, solely in terms of boundary data. 


-/■ 


uj(y) = /[ ti(x) G i: j(x,y) - Fj[ j(x,y) u^x)] dS(x) + / f^(x) G ± j(x,y)dV(x) 


V 


(7) 


It then becomes possible to carry out a limiting operation in which the 
interior point y approaches a specific point x Q on the surface of the body. The 
functions . and Fj. both become singular, and the integral involving F^. 
must be evaluated as the sum of a principal value integral and a separately 


calculated jump term. After the required manipulations have been carried out, 
boundary integral (constraint) equation 


[1 - j(x)] Ui(x 0 ) 






(x) G i j(x,x Q ) - F^Cx.Xq) u^ (x) ] dS(x) 


S 

+ / £ i (x) GijCx.Xo) dV(x) 
V 


the 


( 8 ) 


152 



is obtained. The tensor (X^a is due to the jump term in the integral of F-jj. At 
any point at which the surface is smooth, = 1/2 5 At other points its 

behavior is more complex. 

The basic content of this equation is that a well posed three-dimensional 
elasticity problem can be solved solely in terms of surface displacements and 
tractions. The equations represent constraints on the unknown boundary data which 
ensure that they will be in equilibrium with the prescribed boundary data. 


ELASTIC NUMERICAL IMPLEMENTATION 


While the analytical form of the boundary integral equation is of considerable 
theoretical interest, there are almost no geometry and loading combinations of any 
practical interest for which an exact solution exists. Classical forms of 
approximation, such as series solutions, are cumbersome and, often, of questionable 
accuracy. Practical exploitation of the equation thus requires a technique for 
reducing it to a linear algebraic system. 

Since the boundary integral equation is written in terms of the displacements 
and tractions on the surface of the body, it is necessary to define a convenient 
geometric approximation to the true surface and to define approximations to the 
displacement and traction variation over this surface. This allows the evaluation 
of the required integrals in terms of the algebraic parameters defining the 
displacement and traction variation and, thus, reduction of the overall system to a 
set of linear algebraic equations. If nonzero body forces are involved, then an 
appropriate interior approximation may also be required. It need not be related to, 
or consistent with the surface approximation. The solution is still carried out 
solely in terms of boundary data. 

Two major approaches have been taken to the three-dimensional numerical 
implementation : 

1. The surface is modeled as an assembly of plane triangular patches. 
Displacement and traction are taken to be either constant (ref. 3) or 
linearly varying (ref. 4) over each such patch. Using these approximations 
it is possible to evaluate all required integrals over each surface patch 
in closed form. This leads to an efficient and highly accurate calculation 
of the coefficients in the resulting algebraic system. However, the 
representation of complex geometry requires the use of a large number of 
elements, and the limitation to linear variation of surface data tends to 
reduce accuracy in cases where significant bending or large surface 
element aspect ratios are involved. 

2. The surface is modeled using quadratic isoparametric (fig. 1) surface 
patches (either quadrilateral or triangular). These are derived from the 
use of similar elements in finite element analysis. Surface data is then 
modeled, over each element, using the linear, quadratic or cubic 
isoparametric shape functions (ref. 5). Experience with this form of 
representation has shown that complex geometries can be modeled using a 
reasonable number of elements. The use of quadratic surface data variation 
yields a significant improvement in accuracy, so that a geometrically 
satisfactory model normally possesses sufficient refinement to give 
sufficient definition of stresses and strains. The use of cubic variation 
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has not been found necessary in modeling gas turbine engine components. 

The use of these more general approximations requires the (at least 
partial) use of numerical integration in evaluating the coefficients of 
the linear equation system. 

Using either form of approximation, it is possible to develop algebraic 
representations for interior quantities. 

While both the boundary element and finite element methods are numerical 
stress analysis techniques, significant differences exist between them. In 
particular: 

1- In the finite element method displacements are approximated throughout the 
body. Stress and strain are always derived quantities. In the boundary 
element method only the surface variations of both displacement and 
traction are modeled. Interior results are calculated only at points where 
they are desired. These interior quantities satisfy the exact equations of 
elasticity. 

2. A single finite element is a structure. That is, an elasticity problem can 
be posed and solved for a single element. A single patch used in boundary 
element method modeling is not a structure. It requires an assembly of 
patches enclosing a volume before an elasticity problem can be solved. 

3. The structure of a finite element system matrix and a boundary element 
system matrix differ. A finite element matrix is typically large (since 
the entire volume must be modeled) but is banded and symmetric. The 
boundary element method matrix is significantly smaller but is typically 
full (because of the surface integrals) and nonsymmetric • 

TURBINE RIM COOLING HOLE 


The boundary element method has been applied (using both the forms of 
approximation discussed above) to the calculation of stresses «and strains at 
turbine cooling hole exits, which are possible life limiting locations for turbine 
disks (ref. 6). A test program was carried out, using a large scale model (fig. 2), 
to determine the strain variation around the exit of a radially oriented hole 
intended to feed cooling air to the blade (fig. 3). The program was designed to 
study the effect of different strategies for blending the sharp corner of the hole 
exit . 


An initial attempt was made to analyze the structure using a three-dimensional 
finite element code. It did not prove possible to resolve the strain concentration 
at the hole exit at a reasonable cost (for either model generation or computer 
time). In order to resolve the local behavior a boundary element analysis, using 
plane triangular elements (the BINTEQ code), was carried out. The quarter section 
shown in figure 3 was modeled, as shown in figure 4, using 436 surface patches. 
Boundary conditions were taken from the finite element analysis, based on the 
assumption that the overall load transfer in the finite element analysis was 
correct. Results from the analysis showed good agreement with the strain gage data 
from the model tests (fig. 5). 
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The boundary element method analysis was later redone using isoparametric 
elements (the BASQUE code). The model (fig. 6) required only 97 patches, a 
reduction factor of 4.5. This reduction is directly related to the use of curved 
elements, particularly in the doubly curved region at the hole exit. The results of 
this analysis showed even better agreement with the data than did those of the 
original analysis using plane triangles (fig. 7). The cost was significantly less, 
although computing system changes make an exact comparison impossible. 

The results of these analyses suggest that the boundary element method is 
particularly effective in resolving rapid stress and strain variations, especially 
when accurate geometrical modeling is also required. This observation is consistent 
with results obtained in other analyses as well (ref. 6, 7). It is not yet clear 
whether the boundary element method will be competitive with finite elements for 
the solution of general three-dimensional load transfer problems not involving high 
gradients. In general, the methods should be viewed as complementary, not as 
competitors, of which only one will eventually be used. 


INHOMOGENEOUS AND INELASTIC ANALYSIS 


It was noted above that the analysis of hot section components typically 
involves centrifugal and nonsteady thermal loads, as well as both inhomogeneous and 
inelastic material response. Both the centrifugal and thermal loads can be 
incorporated directly as body forces. For the centrifugal loads (where R is the 
distance from the axis of rotation) 


fi cc pu 2 r 


(9) 


while for the thermal loads 


f 


i 


- aE 

(1 - 2 v ) 


0,1 


( 10 ) 


It should be noted that these representations are not unique. The use of the 
divergence theorem allows transformations of the body force terms to be carried 
out. While the resulting representations are analytically equivalent the choice 
made can have significant numerical effects. 

The introduction of inhomogeneous and inelastic response is more complex. The 
homogeneous, elastic boundary integral equation can be written, in operator 
notation, as 


L(v) = F 


( 11 ) 


155 



where L is a linear operator, v is the unknown surface data and F is derived from 
all known loads and surface data. Other linear operators exist for the calculation 
of any desired interior results: 


u = I U (V) + F u 
e = i e (v) + F e 
O = I a (V) + f ct 


( 12 ) 


The stresses and strains obtained using the operators I e and I a will satisfy 
Hooke's Law using the assumed homogeneous elastic modulus. If either inhomogeneity 
or inelastic response occurs, then the interior stresses and strains will not 
satisfy the true stress strain relation. This imbalance can be represented in terms 
of modified body forces and surface tractions or as volume distributions of initial 
stress or initial strain. If, for example, the body force representation is chosen, 
then the imbalance can be cycled back to the boundary integral equation 


[ i - «i 3 J 


U J - f‘l 


G ij - F i; 


1 

u i 


] dS + 


/ fi 


fi 


1 G ij 


dV 


( 13 ) 


where is derived from the interior imbalance evaluation. The new boundary 
solution u£, will produce interior results which reduce the interior 
imbalance due to the inhomogeneous and inelastic behavior. In matrix form the 
system equation becomes 


[T] u 1 + [U] t 1 = £_ + f 1 


( 14 ) 


It is particularly important to note that the kernel functions and the 
matrices [T] and [U] remain unchanged during the iteration process. Thus the 
calculation of these matrices and the decomposition of the system matrices need be 
done only once for the entire analysis. The overall solution sequence, written in 
operator form is then 

calculate L 
calculate N 
calculate L“l 

calculate F 
f° = 0 

r —v 1 = L -1 (F 
fi+i = n (V 1 ) 

t = t + At 


single elastic calculations 


( 15 ) 


+ f 


f) l. inelastic/inhomogeneous 

/ iteration 


time/load 

stepping 
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The sequence presented above is valid for a single time or load increment. The 
quantities u, t and F must, for an inelastic analysis, be regarded as either 
increments or rates depending on the type of inelastic material model used. In 
either case a number of time steps or load increments will be required for the 
complete solution. The key point to be noted is that the surface and volume 
problems have been uncoupled. In a sense the overall structural analysis problem is 
solved on the surface while the complex material response is evaluated in the 
interior. It remains unnecessary to have compatibility between the surface and 
volume discretizations. 

Significant progress has been made in the application of these techniques in 
two-dimensional analysis (ref. 8, 9). It is to be expected that, with careful 
attention to the details of the numerical implementation, the inelastic 
three-dimensional boundary element method will prove to be highly effective for 
many problems of practical importance. 
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Figure 4 Piecewise Linear Three-Dimensional Boundary Integral Equation Models for 
Analysis of Turbine Disk Rim Slot-Cooling Hole Intersection 
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Figure 5 Correlation of Test and Boundary Integral Equation Results-Second 
Intersection Geometry, Planes 2A and 4A 
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Figure 6 


Figure 7 
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Quadratic Isoparametric Boundary Integral Equation Map for Analysis of 
Turbine Disk Rim Slot-Cooling Hole Intersection 
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Improvement Due to Higher Order Boundary Integral Equation Modeling 
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SUMMARY 


Christos C. Chamis 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135 

Significant results, conclusions, and issues addressed at the Lewis/ 
Industry/University Workshop on Nonlinear Structural Analysis are summarized 
below by discipline: New Concepts/Formulations, Algorithms/Convergence, and 

Inelastic Analysis and Interactive elements. 

New concepts/formulations . - (1) The slave finite element approach to 
nonlinear analysis is viable. The large number of degrees of freedom (DOF) 
associated with this approach, however, would require innovative computational 
algorithms to obtain computational effectiveness compared with traditional 
approaches. The ease in changing increment steps, the increased detailed 
information available within the subinterval, and accuracy demonstrated to 
date for elastic problems are promising. (2) A new variational principle has 
been developed. This principle is based on the combined Hellinger-Reissner 
and Hu-Washizu principles with additional internal displacement DOF. Hybrid 
stress elements based on this new principle have invariant properties and are 
less sensitive to geometric distortions. This new variational principle is 
considered to represent a significant advance in finite element development. 

(3) The three-dimensional degenerated elements offer modeling and computa- 
tional efficiency for the analysis of general shells with layered structure, 
undergoing large deformations (in-plane, bending, and shear) and transient 
motions. As expected, these elements take more computational time for special 
shell problems compared with finite elements based on deformable shell 
theories. (4) A nine-node Lagrange shell finite element proved effective in 
the buckling, postbuckling, and elastic-plastic analyses of shells. This ele- 
ment did not exhibit "shear locking" even in very large aspect ratio shells 
(up to 10 8 ). Application of the element to other nonlinear problems such as 
dynamic, creep buckling, and ratcheting requires additional research. The 
magnitude of the effort required to develop the element formulation for these 
other nonlinear problems may be substantially reduced using symbolic language 
manipulators. 

Algorithms/convergence. - (1) The new computational algorithm developed 
for self-adaptive solution strategies proved to be very effective in the solu- 
tion of two highly nonlinear problems (snap-through load-deformation history 
of a spherical cap and a box-truss). This new algorithm is stable, eliminates 
the need for global inverse, and handles problems with turning points and 
general material properties. Formal establishment of convergence and stabil- 
ity characteristics of the algorithm requires additional research. (2) 
Element-by-element (EBE) approximate factorization techniques were compared to 
test problems for elastic and elastic-perfectly plastic cantilever beams with 
and without a circular hole. The EBE solution method converged with seven 
iterations (with the hole) and four iterations (without the hole). EBE is a 
linear equation-solving procedure. Direct use to develop EBE time stepping- 
solution algorithms requires additonal research. (3) Finite element genera- 
tors using symbolic manipulator languages (VAXIMA) are viable. Initial 
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results show that the lines of program required to generate the strain- 
displacement and stiffness element matrices can be reduced substantially. (4) 
The stability convergence of underintegrated finite element approximations can 
be improved by simple post-processing operations. These operations can be 
described by an exact formula in some cases. Additional research is needed 
for more general cases. 

Inelastic analysis and interactive elements . - (1) Elastic-plastic frac- 
ture under arbitrary loading and material behavior, characterized by realistic 
crack-growth rates, can be described by two integral equations. A new 
path-independent integral has been derived to describe dynamic crack propaga- 
tion. A set of equations/conditions has been developed which generalizes the 
constitutive relations for inelasticity and finite deformations. Variational 
principles were used to derive a hybrid finite element for describing dynamic 
crack propagation. (2) Interactive finite elements can be used to describe 
the dynamic interaction of rotating squeeze film bearings, shafts, and sup- 
ports. This type of analysis is the only one available to date which 
describes the dynamic interaction in an integrated manner and which can be 
used as a "plug-in" in general purpose finite element codes. The interactive 
finite element approach is presently being confirmed experimentally using a 
two-bearing, two-damper test rig. (3) An upward-integrated and "top-down 
traced" computerized capability has been developed for the nonlinear analysis 
of high-temperature composites applicable to turbine blades and vanes. This 
theory accounts for the complex thermoviscoplastic behavior of composites in 
turbine environments and includes effects of the fabrication process. Initial 
results indicate that the fabrication process induces a substantial thermo- 
viscoplastic stress state in the composite. (4) Boundary element methods have 
proven computationally effective for two-dimensional linear problems. Exten- 
sions to three-dimensional, nonlinear and time-dependent problems is the 
subject of current focused research. It is anticipated that boundary elements 
will provide a complementary analysis methodology to finite element methods. 

Collectively, the summary points out the multidisciplinary approach 
necessary to develop methods for complex nonlinear structural analysis 
problems. This multidisciplinary approach encompasses new formulations and 
concepts, invention of more general variational principles, finite element 
generators, development of interactive elements, innovative solution algo- 
rithms, adaptive solution strategies, and formal criteria for a priori assess- 
ing solution stability and convergence. The contributors to these proceedings 
describe what can be done, future potential, and direction in each discipline. 

The goal of the highly focused research described in these proceedings is 
to develop methodologies which lead to designs of more cost-effective and 
reliable engine structures and components. It is hoped that, in addition, 
these proceedings will provide a valuable source of information of on-going 
research and future direction in nonlinear structural analysis in general. 
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INDEX 


Acceleration, secant, 48 
Adnissible, 3 
Admissible, 3 
Ahmad, 45, 49 
Algorithm, 

Hierarchial, 58 

Hybrid-stress finite element, 111 
Incremental Newton-Raphson, 55 
Newmark, 72 

Newton-Raphson, 46, 49 
Parabolic regularization, 66, 67 
Analysis 

Boundary element method, 154 
Global 132 
Postbuckling, 46 
Thin plate, 47 
Thin shell , 47 
Top-down-structured 132 
Upward-integrated 132 
Application levels, hierarchial, 57 
Arc length, constrained, 46, 49 
Aspect ratio. Large, iii 
Attributes, self-adaptive, 56, 57 

Beam, cantilever, 71 

Bearing dampers, squeeze-film, 121 

Bearings 

Fluid-film journal, 120 
Rolling element, 120 
Behavior 

Inelastic, 156 
Inhomogeneous, 156 
Pre-postbuckling, 60 
Rate sensitive, 106 
Bending 

Modes, inextensional , 20 
Static, 35 
Bilinear, 97 

Blade/vane components, high-pressure turbine 
(HPT), 131, 132 

Boundary element, 154, 155, 162 
Boundary element method, 149 
Boundary integral, iii 

Cap 

Sandwich spherical, 48 
Spherical, 36, 60 
Closedness criterion, 56 
Composites 

Constituent materials, 132 - 134 
Fiber reinforced, 31 
High-temperature materials, 131, 162 


Composites (cont.) 

Inhomogeneous TFRS, 134 
Intralaminar 
Hicrostresses, 133 
Non-uniformity, 134 
Lamina, 133 

Mechanical strength, 133 
Mechanics, 133 
Micromechanics, 133, 134 
Refractory-fiber metal-matrix, 131 
Thermal strength, 133 
Unidirectional TFRS, 134, 135 
see Tungsten-fiber-reinforced superalloys 
Computer code 
ADINA, 122 

COBSTRAN, 132, 133, 134 
FORTRAN, 85, 86 
FORTRAN 77, 88, 89 
Interactive prompting, 87 
LISP, 87, 89, 90 
MACSYMA, iii, 85, 86 
Menu driven, 87 
NASTRAN, 120, 121, 134 
TETRA, 120 

User-level programs, 90 
VAXIMA, 86, 87, 89 

Computer operating system, UNIX, 86, 88 
Computers, VAX-11/780, 86, 91 
Concave shape, 49 
Configuration, multi-shaft, 120 
Conservation laws, 108 
Constitutive laws, 3 
Constitutive relation, 105 
Constraint 
Circular, 56 
Hyperbolic, 56 
Surface of, 56 
Contour line, 74 
Convergence, 162 
Criteria, iii 
Rates, super 1 inear, 56 
Stable, 57 
Convex shape, 49 
Coordinates, curvilinear, 33 
Cowper, 22 
Crack propagation 
Dynamic, 105, 162 
Elastodynamic, 108 
Self-similar, 1052, 108 
See also creep 
Creep 

Behavior, elastic-plastic, 55, 56 
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Creep (cont.) 

Crack-growth, 106 
Damage process, 106 
Steady-state 106 
Cycle, limit, 121 

Dampers, squeeze-film, 120, 122, 123 
Dawe, 22 

Deflection, large, 48 
Deformation 

Elastic-plastic, 48 
Inelastic, 112 
Large, 46, 47 
Theory of, 105 
Determinant ratio, 49 
Derivative 
Spatial, 3 
Temporal, 3 
Discretizations, 65 
Non-prismatic, 10 
Semi -prismatic, 10 
Displacements 
Constrained, 48 
Internal, 20 

Statically condensed, 18 
Distortions, geometric, 20 
Dynamics, nonlinear, 121 

Eight-node, 20 

Elastic-plastic material, Bilinear, 48 
Elastic-plastic, 46 
Elasticity 
Finite, 106 
Infinitesimal, 106 
Element-by-element, iii 
Elements 

Bilinear quadrilateral, 72 
Boundary, 151 
Cylindrical shell, 22 
Degenerated, 161 
Higher order, 87 
Hybrid, 17 
Hybrid semi-Loof, 21 
Implicit-explicit, 71 
Invariant, 17, 20 
Isoparametric, 155 
Lagrange, 46 
Nine-node, 47 

Quadralateral isoparametric, 60 
Semi -1 oof, 17,20-22 
Thin shell, 20, 21 

Three-dimensional degenerated, 31, 32, 40 
Time-adaptive implicit-explicit partition 
of, 71 

Two-dimensional, 40 


Elliptic, 56 

Energy, complementary, iii, 112 
Equation 

Boundary integral, 153 
Euler-Lagrange, 112 
Micrcmechanics, 132 
Euclidean norm, 58 
Euler-Lagrange, 112 
Factorization, 76 

Alternating direction method, 65 
Element-by element, 65, 66, 70, 71, 161 
Fractional steps method, 65 
Methods of discretization, 65, 66 
Splitting-up method, 65 
Weak approximation method, 65 
Finite elements, 33, 85, 86, 111 
Analysis, 45, 55 
Conforming, approximation, 97 
Generators of, 161 
Global, 134 
Hybrid stress, iii 
Interactive, iii, 162 
Nonlinear simulation, 55 
Slave, iii, 161 

Flow theory, rate-independent, 106 
Force excitation, periodic, 121 
Formulation 
Galerkin, 57 
Rayleigh-Ritz, 57 
Virtual work, 57 
Fracture 

Dynamic, iii 
Mechanics, 105, 107 
Elastic-plastic, 107, 162 
Frequencies, natural, 121 
Function, bilinear, 7, 8 
Functional, 112 

Global inversion, 57 
Global updating, 57 
Green-Lagrange, 32, 33 
Growth, self-similar, 105 

Hamilton's law, 1, 2 
Hardening, anisotropic, 109 
Heat conduction, nonlinear, 56 
Heirarchy, 132 
Hellinger-Reissner, 17, 18 
Higher order, 46 
Hinge, plastic, 48 
Hu-Wash izu, 17 
Hybrid elements, 17 
Hypo-elasticity, 109 

Incompressibility, 95 
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Inelastic analysis, iii, 149 
Instabilities, hourglass, 97 
Instantaneous, 122 
Integral 

Path-independent, 105, 106, 108, 109, 162 
Principal value, 152 
Integration 
Reduced, 45 
direct, Newmark, 35 
Interphase, 132, 135 
Invariant, 18 
Iron, 45, 49 
Isoparametric, 18 
Instabilities, numerical, 95 
Iteration, Newton-Raphson, 48 
J integral, 105 
Jaumann integral, 115 
Jump term, 152, 153 

Labyrinths, 121 
Lagrange, 46 
Lagrange element 
Higher order, 47 
Nine-node, 48 
Lagrangian 

Kinematical, 75 
Total, 32 

Lamina transverse directions, 133 
Laminate theory, 132, 134 
Large global arrays, factorizing, 65 
Limit cycle, 121 
Load, temperature, 135 
Loading/unloading histories, 106 

Manipulation, algebraic, 85 

Marguerre, 22 

Material 

Elastic, 48 

Elastic-plastic, 48, 109 
Elastic-perfectly plastic, 48, 72, 73 
Rigid-plastic, 109 
Matrix 

Lumped mass, 72 
Tangent damping, 77 
Tangent stiffness, 77 
Membrane 
Action, 49 
Quadrilateral, 18 
Stresses, momentless, 22 
Mesh partitions, implicit-explicit, 66 
Micro-mechanics, 109 
Mindlin, 46 
Modes 

Checkerboard, 95 
Chickenwire, 95 


Modes (cont.) 

Free-free undamped normal, 121 
Hourglass, 95, 98, 99 
Keystone, 95 
Rigid body, 20 
Morley, 22 

Motion, geometrically nonlinear, 32 
Multi-branch, multi-level system, 120, 122 
Multi -component splitting formulae 
Two-pass, 69, 70 
One-pass, 70 

Multiprocessor computer, 66 

NFAP, 88, 91 
Newmark, 35, 75 
Newton-Raphson, 46 
Nonlinear, 85 
Nonlinearity 
Kinematic, 55 
Material, 55 
Noor, 36 

Normality condition, 109 

Objectivity, 113, 105, 114 

Panel, circular cylindrical, 35 
Parallelizable, 70 

Partitions, implicit-explicit mesh, 77 

Patch test, 18 

Patches 

Quadratic isoparametric surface, 153 
Triangular, 153 
Path-independent, iii 
Patterns, spurious, 95 
Phenomenon, transient, 121 
Piola-Kirchoff . 32, 58 
Plastic deformation, non-proportional, 105 
Plastic zone, 49 
Plasticity theory 

Linear kinematic hardening, 110 
Isotropic-hardening, 109 
Plates, twisted, 36 
Polynomials 
Bilinear, 18 
Interpolation, 88 
Quadratic, 46 
Positive definite, 71, 77 
Positive semi -definite, 77 
Postbuckling, 48, 49 
Algorithm, 48 

Analysis, artificial spring method, 46 
Prebuckling, 49 
Pre-postbuckling, elastic, 56 
Pressure, 48 
Follower, 48 
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Pressure (cont.) 

Constant-direction, 48 
Principle, variational, iii, 18, 161 
Prismatic, 2 

Propagating singular element, 108 

Quadrature 
Gaussian, 98 
Numerical, 95 

Quadrilateral planform, 21 

Ramberg-Osgood, 7 
Ramm, 46, 50 
Rank, 18 
Rate 

Objective, 109, 110 
Strain, 111 
Stress, 110, 112 
Zaremba-Jaumann-Nol 1 , 111 
Rayleigh-Ritz, 58 

Relations, macroscopic constitutive, 109 

Resonance, subharmonic, 121 

Response, transient, 35, 36 

Riks, 46, 50 

Rolling contact, 56 

Rotor-stator 

As interactive components, 120 
Clearances, 121 
Interactive forces, 122 
Rub/impact, 122 

Safety zones, 56 
Seals, 122 

Secant accelerated, 46 
Self-adaptive, iii 
Semi-discrete, 75 
SemiLoof, see Elements 
Sensitivity studies, 132 
Serendipity, 46 
Shape function, 87 
Shear 

Locking, 45, 47 
Finite, 111 
Shell 

Angle-ply, 36, 40 
Cross-ply, 36, 40 
Cylindrical, 36, 40 
Degenerated, 45, 47 
Elements, 47 

Laminated anisotropic, 31 
Large aspect ratio, 161 
Shear-deformed, iii 
Spherical, 36 
Spherical shallow, 48 
Theory, 22 
Thin rotating, 120 


Shell theory 
Cylindrical, 22 
Deep shallow, 22 
Shallow, 22 

Solid element, degenerated, 22 
Solution 
Curve, 56 
Self-adaptive, 161 
Somigliana's identity, 152 
Space-time, 3 

Splitting, two-component, 69 
Stationarity, 112 

Stiffness, incremental tangent, 46 
Strain 

Green-Lagrangian, 33 
Momentless membrane, 20 
Rate, 111 
Stress 

Dependency, 135 
Hybrid, 19 
Incremental, 33 

Large deformation inelastic analysis, 105 
Tensor, Piola-Kirchhoff, 32, 58 
Stress-strain field. Crack-tip, 108 
Stricklin, 36 
Structural analysis, 134 
Structure, box truss, 60 
Subdomain model, 70 
Subgroups, Non-conti guous, 70 
Substructures, 70 
Superelements, 70 
Swept panel, 20 
Symbolic, 85, 86 
Expressions, 89 
Processing, 87 

TFRS, see Tungsten-fiber-reinforced superalloys 

Tangent stiffness, singularity of, 46 

Taylor series, 58 

Temperature dependency, 135 

Temporal , 3 

Test rig, 123 

Theory 

Deformation, 105 
Four-dimensional unified, 10 
Laminate, 132 

Thermoviscoplastic (TVP), nonlinear, iii, 132, 
133, 135 

Thin plate analysis, 47 
Thin shell analysis, 47 
Top-down traced, 162 
Transient, geometrically nonlinear, 31 
Triangular elements, 21, 22 
Tungsten-fiber-reinforced superalloys (TFRS), 
131, 132, 134 
Blades/ vanes, 133 
Fiber degradation of, 132 
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Tungsten-fiber-reinforced superalloys (cont.) Vibration, natural, 35, 36 

Fiber-matrix interface of, 132 Viscoplasticity, 106 

Nonlinear analysis of, 132 

Turning points, 56 Work theorem, reciprocal, 152 


Underintegration, 95 Zienkiewicz, 45, 49 

Bilinear, 98 Zone, degraded, 135 

Unstable, 48 
Upward-integrated, 162 
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